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ABSTRACT 

To  aid  in  the  modeling  of  a  steady  state  spin,  the 
equations  of  motion  of  an  airplane  are  formulated  in  a 
cylindrical  coordinate  reference  frame.   The  derivation 
of  the  equations  is  presented  and  the  resulting  equa- 
tions are  simplified  for  the  equilibrium  spin  condition. 
These  simplified  equations  are  used  in  an  unconstrained 
computer  parameter  optimization  technique  that  algebrai- 
cally solves  the  differential  equations  for  the  equilib- 
rium state.    The  results  of  the  computer  work  are 
presented  and  compared  with  previous  prediction  schemes, 
The  potential  of  the  method  is  demonstrated  by  applica- 
tion to  a  study  of  the  effects  of  density  variation. 
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LIST  OF  SYMBOLS 

The  definitions  of  the  symbols  used  in  this  paper  are  as 

follows : 

b      Wing  span,-  ft 

c      Mean  aerodynamic  chord,  ft 

F_       Lagrange  generalized  force,  lbs   or 
4r      ft-lbs 

F  , F  , F      Aerodynamic  forces  along  the  body  axes, 
lbs 

2 

g     Acceleration  of  gravity,  ft/sec 

I  ,1  ,1  ,1        Body  axes  moments  and  products  of 

inertia,  slugs— ft 

M      Rolling  moment  acting  about  the  X 
body  axis,  ft— lbs 

M       Pitching  moment  acting  about  the  Y 
body  axis,  ft— lbs 

M       Yawing  moment  acting  about  the  Z 
body  axis,  ft— lbs 

m     Airplane  mass,  slugs 

p,q,r      Angular  rates  of  rotation  about  the 
X,  Y  and  Z  body  axes,  rad/sec 

2 
q^      Dynamic  pressure,  1/2  pV 

q       Lagrange  generalized  coordinate 

R      Radius  of  helical  path  of  aircraft,  ft 
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S     Wing  area,  ft 

T  Lagrange  kinetic  energy,  ft— lbs 

V  Inertial  velocity  of  the  C.G.,  ft/sec 

v  Lagrange  potential  energy,  ft— lbs 

X,Y,Z  Aircraft  body  axes 


X  ,Y  ,Z,      Cartesian  coordinate  reference  axes 
that  are  attached  to  the  cylindrical 
radius  vector. 

Z       Aircraft  altitude,  ft 
o 

a  Angle  of  attack,  degrees 

a.  .  Direction  cosine,  dimensionless 

3  Angle  of  sideslip,  degrees 

Y  Cylindrical  orientation  angle,  radians 

6       Aileron  deflection  (positive  when  right 
aileron  trailing  edge  is  down) ,  degrees 

6^     Elevator  deflection  (positive  when 

elevator  trailing  edge  is  down) ,  degrees 
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6       Rudder  deflection  (positive  when  rudder 
trailing  edge  is  left  when  viewed  from 
above) ,  degrees 
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p      Atmospheric  density,  slugs/ft 

ip,6,<}>      Euler  orientation  angles 
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I.   INTRODUCTION 

A.   BACKGROUND 

As  a  result  of  post— stall/spin  accidents  a  significant 
number  of  lost  fighter  aircraft  and  pilot  fatalities  have 
been  recorded  over  the  last  decade.   Investigations  have 
shown  that  these  aircraft  exhibit  poor  spin  characteristics 
and  that  recovery  from  a  fully— developed  spin  is  usually 
difficult  or  impossible.   For  these  reasons  the  spin  is  no 
longer  considered  a  tactical  maneuver  and  is  regarded  by 
pilots  as  an  undesirable  and  potentially  dangerous  flight 
condition. 

Although  the  flight  conditions  leading  to  a  spin  can 
generally  be  avoided,  the  combat  pilot  often  finds  the  need 
to  utilize  the  edge  of  the  flight  envelope  in  a  combat 
situation.   More  often  than  not,  it  is  only  in  this  boundary 
area  that  the  combat  superiority  of  one  aircraft  over  an- 
other is  realized.   If  the  aircraft  demonstrates  unsatis- 
factory stall  and  spin  characteristics,  the  pilot  is 
disinclined  or  restricted  from  utilizing  this  outermost 
region  of  flight  to  his  best  advantage.   As  a  result,  the 
overall  weapon  system  effectiveness  and  margin  of  superi- 
ority of  the  aircraft  is  lost. 

In  light  of  the  poor  stall  and  spin  characteristics 
of  current  fighter  configurations  a  more  detailed  consid- 
eration must  be  given  to  these  characteristics  during  the 
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early  design  stages  of  future  fighter  aircraft.   In  recog- 
nition of  this  fact  this  research  was  undertaken  to  provide 
a  more  accurate  analytic  tool  with  which  the  designer  can 
investigate  the  spin  characteristics  of  future  aircraft  and 
to  aid  in  understanding  how  various  factors  in  design  may 
affect  these  spin  motions.   It  was  not  proposed  that  this 
effort  solve  the  spin  problem  itself,  but  it  was  hoped  that 
it  would  prove  to  be  a  useful  approach  and  that  it  would 
provide  information  complementary  to  that  obtained  experi- 
mentally. 

B.   SCOPE 

Traditionally,  analytic  spin  prediction  has  involved 
either  an  approximated  solution  (by  requiring  that  only 
part  of  the  equilibrium  conditions  be  satisfied)  or  by  com- 
puter generation  of  time— histories  of  characteristic  para- 
meters (by  integrating  the  equations  of  motion  forward  in 
time  until  the  average  values  of  the  variables  are  approxi- 
mately constant) .   More  recently  the  idea  of  mathematically 
casting  the  equations  of  motion  into  a  form  specifically 
intended  for  spin  prediction  has  shown  some  promise.   Of 
particular  note  is  the  work  carried  on  by  Buehler  (Ref.  1) 
and  Champoux  (Ref.  2).   These  earlier  works  provided  much 
of  the  initial  groundwork  for  this  research. 

In  developing  the  analytic  approach  to  the  spin  problem 
a  set  of  objectives  was  established  to  subdivide  the  re- 
search problem  and  also  to  serve  as  guidelines  to  ensure 
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that  the  work  proceeded  in  the  proper  direction  as  effec- 
tively as  possible.   These  objectives  were: 

1.  To  suitably  model  the  aircraft  in  a  spin. 

2.  To  develop  the  spin  equations  of  motion. 

3.  To  establish  a  suitable  technique  of  solving  the 
equations  of  motion  for  the  equilibrium  spin  con- 
dition . 

4.  To  develop  a  computer  program  that  will  predict 
aircraft  spin  characteristics.   In  addition,  the 
program  should  be  flexible  enough  to  aid  in  the 
understanding  of  how  various  parameters  effect  the 
developed  spin  motions. 

5.  To  demonstrate  the  usefulness  of  the  method  and  to 
verify  the  results  with  other  prediction  schemes. 
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II.   DEVELOPMENT  OF  THE  AIRCRAFT  EQUATIONS  OF  MOTION 

A.   COORDINATE  SYSTEM 

To  an  observer  on  earth,  the  aircraft's  center  of  grav- 
ity in  a  steady  spin  appears  to  trace  a  descending  constant 
radius  helical  path  about  an  imaginary  vertical  central 
axis.   By  assuming  that  the  rotation  of  the  earth  is  negli- 
gible, an  inertial  reference  frame  can  be  established  that 
is  fixed  to  the  earth's  surface  at  the  point  of  intersection 
of  this  imaginary  spin  axis. 

The  most  convenient  method  of  describing  this  motion  is 
through  the  use  of  cylindrical  coordinates.   In  this  formu- 
lation, the  central  axis  of  the  cylinder   is  superimposed 
on  the  imaginary  aircraft  spin  axis  as  shown  in  Figure  1. 
The  advantage  of  this  choice  for  an  inertial  coordinate  sys- 
tem is  that  the  center  of  gravity  of  the  aircraft  is  con- 
veniently located  in  terms  of  an  altitude  coordinate  (Z  ) , 
a  radial  coordinate  (R)  and  an  angular  position  coordinate 

(Y). 

The  body  axes,  as  shown  in  Figure  2,  were  chosen  for 
the  formulation  of  the  equations  of  motion.   This  allows 
the  direct  use  of  standard  wind  tunnel  data  which  was  avail- 
able- for  this  research  effort.   Symmetry  of  the  aircraft 
about  the  X  and  Z  axes  was  assumed  so  that  all  cross— product 

of  inertia  terms,  with  the  exception  of  I   ,  were  identi— 

XZ 

cally  zero. 
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Figure  1.   Inertial  Cylindrical  Coordinate 
System. 


X 


.— e 


Figure  2.   Aircraft  Axes 
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The  orientation  of  the  aircraft  is  described  by  intro- 
ducing another  coordinate  system  and  a  set  of  three  rota- 
tion angles.   A  fixed  cartesian  coordinate  reference  frame 
(X, ,Y.,Z ,)  is  positioned  at  the  center  of  gravity  of  the 
aircraft  and  aligned  with  respect  to  the  inertial  cylindri- 
cal frame  so  that  the  Z,  axis  is  parallel  with  the  vertical 
central  axis  and  the  Y,  axis  coincides  with  the  — R  position 
vector.   Although  this  cartesian  coordinate  reference  frame 
remains  fixed  in  orientation  with  respect  to  the  inertial 
cylindrical  frame,  it  is  still  free  to  move  in  space  as  Z  , 
R  and  y    take  on  different  values.   The  orientation  of  the 
body  axes  relative  to  the  cylindrical  frame  of  reference  is 
given  by  the  aircraft  Euler  angles  (ip,8,(f)),  which  define  an 
ordered  series  of  three  consecutive  rotations  as  shown  in 
Figure  3.   The  three  position  coordinates  (Z  ,R,y),  together 
with  the  three  orientation  angles  {ty,Q,$),    allow  a  complete 
description  of  the  motion  of  the  aircraft  in  space  and  will 
serve  as  the  generalized  coordinates  in  the  Lagrangian 
development  of  the  aircraft  equations  of  motion. 

B.   AXES  TRANSFORMATION  AND  AIRCRAFT  ANGULAR  RATES 

Before  proceeding  with  the  development  of  the  equations 
of  motion,  certain  algebraic  relations  are  needed  that  will 
transform  the  actual  forces,  moments,  velocities  and  angular 
rates  of  the  aircraft,  about  its  body  axes,  to  the  cylindri- 
cal reference  frame.   These  transformations  are  nothing  more 
than  vector  projections  from  one  system  of  axes  onto  another. 
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ihertiaL 


Figure  3.   Spin  Model. 

Consider  that  the  body  axes  are  initially  coincident 
with  the  inertial  frame  (X, ,Y. ,Z.)  and  then  are  rotated 

about  Z,  by  an  angle  \p    to  form  an  intermediate  reference 

p 

frame  (X„,Y„,Z„),  as  shown  in  Figure  4 

Since  both  of  these  orthogonal  co- 
ordinate systems  form  a  basis  for 
three— space,  any  point  (P)  can  be 
described  by  a  vector  (R)  from  the 
origin  to  that  point  in  terms  of 
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either  system. 

Mathematically,  this  can  be  expressed 


r  =  X^  +  Y1j1  +  Z^  =  X2i2  +  Y2j2  +  Z2k2 


(1) 


but, 


12"  r 


X2  =  ¥2"  *!  +  YlV  jl  +'  Zli2'  kl 


V  r 


=  Y 


2  =  XlV  £1  +  Ylj2*  jl  +  ZlV  kl 


(2 


k2-  r 


z2  -  x1k2-  ix  +  Ylk2.  j  +  z  k2-  ^ 


which  is  written  more  compactly  in  matrix  form  as 


X. 


V1! 


j2-ix 


k2-i 


/\  /V 


VJ1  W 


/\  <A. 


Vji      Vki 


k2.j1    k  .k 


X. 


zl 


(3 


The  three— by— three  matrix  of  unit  vectors  is  called  a  trans- 
formation matrix  and  denoted,  in  this  case,  by  [T-,,]  and 
read  "the  transformation  from  frame  one  to  two."   Carrying 
out  the  indicated  dot  products  for  the  rotation,  the  trans- 
formation matrix  has  the  value: 


[Tvi] 


cos    \p        —sin    \p  0' 

sin    \p  cos    \p  0 

0  0  1 


(4) 


The  second  rotation,  9,  about  the  Y_  axis  carries  sub- 
scripted coordinate  system  two  into  system  three,  as 
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shown  in  Figure  3.   The  resulting  transformation  matrix 
is : 


[T3/2] 


cos  9 

0 

0 

1 

sin  9 

0 

sin  9" 

0 
cos  { 


(5) 


The  final  rotation,  $ ,  about  the  X-.  axis  carries  sub- 
scripted coordinate  system  three  into  the  unsubscripted 
body  axes,  as  shown  in  Figure  3.   This  resulting  trans- 
formation matrix  is: 


[V3] 


1 

0 

0 

0 

cos  $ 

sin 

<J> 

0 

,— sin  <f) 

cos 

<f> 

(6) 


These  individual  transformation  matrices  can  be  multiplied 

together  to  form  a  single  three— by— three  matrix  which  will 

be  denoted  by  [T  ., ] .   Since  each  transformation  matrix 

p/l 

premultiplies  the  vector  arrived  at  in  the  previous  step, 
the  total  transformation  from  the  cartesian  coordinate 
frame  to  the  aircraft  body  axes  is  given  by: 


[V] 


T 


/3 


x 


3/2 


x 


2/1 


(7) 


Carrying  out  the  matrix  multiplication,  the  total  trans- 


formation has  the  following  form  in  which  a. .  denotes  its 


elements : 


[Vi] 


a,  ,  =  cosBcosiJ; 


a„  =  sin<£>sin9 

COSlJ) 

— cos<J>siniJ; 
oc_    =    cos(|)sin9 

COSl|i 

+sinc})simp 


a,    =    cos9siniJj 


a„p=    sincfjsinO 
sinijj 

+  COS(f)COSlp 

cosd>sin9 


a 


32' 


sin^ 
— sincfjcostj; 


a,    =   — sint 


a„    =    sin<{)COs6 


a,.-.=    coscJjcosi 


(8) 


In  relating  the  angular  velocity  rates  of  the  aircraft, 
p,  q  and  r  to  the  inertial  frame,  we  make  use  of  the  trans- 
formations [T  ,..  ]  ,  [T  ,,.]  and  [T-,,,].   This  development 
entails  projecting  the  orientation  angular  rates,  \p ,    6  and 
(J>  onto  the  aircraft  body  axes.   In  addition,  the  inertial 
angular  rate  (y)  must  be  included  to  account  for  the  rota- 
tion of  the  cartesian  coordinate  reference  frame  (X, ,Y  ,  Z,). 

Both  y  and  i>    are  considered  to  be  applied  about  the  Z, 
axis  of  the  cartesian  coordinate  reference  system  and, 
through  the  [T  ,.  1  transformation,  their  projections  are 
reali2ed.   The  Euler  pitch  rate  (9),  about  the  Y„  axis,  in- 
volves only  two  transformations  [T~/?]  and  [T  ,A  ,    for  pro- 
jection on  the  body  axes.   The  Euler  roll  rate  (<p)    about 
the  body  axis  involves  no  transformations.   Mathematically 
expressed: 


(9 


p 

0 

0  ' 

<$> 

q 

L    P/1J 

0 

+    tTp/3!1T3/21       ' 

8 

1      + 

0 

r 

i  ty+y 

. 

0 

0 
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On  substitution  of  Equations  (5),  (6)  and  (8)  for  [!',„], 

[T  ,-J  and  [T  <,  ]  respectively,  in  Equation  (9),  the  exact 
P/3         p/1 

relation  between  the  aircraft  angular  rates  and  the  Euler 

angle  rates  is  given  by: 

•  •  • 

p  =  (j>  —  sin0  (ip+y) 

q   =    6    coscf)    +    cos6sin<{)  (i^+y)'  (10) 

r   =    cos 0 cost})  (tp+y)    —9    sin;}) 

The  time  rates  of  change  of  the  direction  cosines  (a.  .) 
and  the  angular  rates  (p,q  and  r)  are  obtained  by  direct 
differentiation  and  are  presented  in  Appendix  A. 

C.   LAGRANGE  FORMULATION  OF  THE  AIRCRAFT  EQUATIONS 
OF  MOTION 

Assuming  that  the  aircraft  is  a  single  rigid  body,  the 
three  position  coordinates  (Z  ,R,y)  an<3  the  three  orienta- 
tion angles  (^,0,<f>)  can  completely  describe  the  six  degrees 
of  freedom  of  the  vehicle.   These  six  variables  form  the 
independent  set  of  generalized  coordinates  necessary  for 
Lagrange  modeling  of  the  motion.   Let  T  denote  the  kinetic 
energy  of  the  system;  V,  the  system  potential  energy;  q  , 
a  generalized  coordinate;  and  F_    a  generalized  force  due 
to  non— conservative  forces.   Then  the  Lagrange  equations 
for  this  system  can  be  written: 

-r-r( )  —  -^-  +   -7T-      =  Fq   ,   r  =  1,2, ...6          (11) 

dt  . •  do           dq              ^r 

9q  ^r     ^r 
^r 

The  kinetic  energy  of  the  aircraft,  relative  to  the 
inertial  reference  frame,  is  given  by  the  superposition  of 
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the  energy  of  the  center  of  gravity,  plus  inertial  rota- 
tional terms  about  the  center  of  gravity,  and  is  written: 

T  =  |m[(yR)2+  (R)2+  (ZQ)2]  +  ^[IXP2+  IYq2+  Izr2]+  JXZpr  (12) 

It  is  necessary  to  use  the  relations  given  in  Equation  (10) , 
as  Equation  (12)  is  not  compatible  with  the  generalized  co- 
ordinates (Z  ,R,y  ,\p  ,Q  ,<\>)  .       Upon  substitution  of  these  rela- 
tions for  p,  q  and  r  the  compatible  kinetic  energy  equation 
is  : 

T  =  ^m[(YR)2  +  (R)2+  (Zq)2]  +  |lxti  -  sinG(i+y)]2 

1  •     •        2 

+    -^v^9    COScf)    +    sin4>cose  (^+Y)  1  (13) 

1  •     •       2 

+   tt-I    [cos<J>cos6  (iJm-y)-] 

+    I       \§   —  sin0  (iJi+y)  ]  [(^+Y)    costficosG    —    6    sintfi] 

If   a   reference    for    the   potential    energy    is    established 
at    sea    level,    V   is   given   by: 

V  =   mgZ  (14) 

The  partial  derivatives  of  Equation  (13)  with  respect  to 

the  time  rates  of  the  generalized  coordinates  ( —  )  are: 

3qr 

8T  • 

—  =   mZ  (15) 

9Z 
o 

8T 

—  =   mR  (16) 

8R 

—  =   mRY+   IxPa13+    *YW23+    Va33 

+    I       (p   cos^cose      —  r    sinG)  (17) 
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and 


f?  =    ^^13    +    ZYqa23    +    ZZra33    "    XXZP    Sin*  (18) 

dip 

—  =    I    q    cos  (J)    +    I    r    sine})    +    I       (pcos$cos0— rsin0)        (19) 
30  Y  Z  XZ 

§  =  V  +  xxzq  -  (20) 

9<J> 


The  time  rates  of  change  of  Equations  (15)  through  (20) 

represent  -r—  (—7-  )  and  are: 

3q 
^r 

at  <ff  »  ■  m^o  (21) 

o 

ir    (^?)       =   mR  (22) 

dt      3R 


d       3T  •  •  •  • 

gr;(— 7)    =   mR(RY  +  2Ry) +IX  (pa,3+pa  13)+I  v(qa23+qa23)    + 
3y 

I    (ra.'+ra,,)  +1       (— p  (0cosc})sin0  +  $sin({)cos0)  + 


pcos^cos©—  0rcos0  —  r  sin0)  (23) 


It^    =    VP^13+pa13)+Vq^23+qa23)+V^33+**a33)  + 

dip 

•  •  • 

Iv-7  ("~P  (0cost})sin0  +  ({)sin4)cos0  )  +pcoscj)cos0    — 

0    r    cos0    -  r    sin0)  (24) 


^-(—r)    =    Iy(qcos<|>— q<j>sin<J>)+Iz  (rsin<j>+r$cos<f>) 
30 

— IyZ  (p0cosc{)+psin(t))  (25) 

dt (7;}   =  *xp  +  xxzq  (26) 
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The  partial  derivatives  of  Equation  (13)  with  respect 

3T 
to  the  generalized  coordinates  (^—  )  are: 


'r 


9T 

9Z 
o 

=    0 

9T 
3R 

•2 

=   my    R 

9T 
9y 

=    0 

9T 

=    0 

(27) 
(28) 
(29) 
(30) 


9T  * 

__  =  —(tp+y)  (I    pcos0  +  I    qsin9sin$  +  I    rsmGcosij)) 

a  u  Ax  Z 

+1   (— pcosc{)sin0  (i+y)  -  rcos6(i+y))  (31) 

|^  =  qr(Iy-Iz)  +  lxz(-icos<|)  -  sin0cos6  (i+y)  )  (32) 


The  partial  derivatives  of  Equation  (14)  with  respect 

9  V 
to  the  generalized  coordinates  (tc—  )  are: 
^  9q 

^r 

ff-  =  mg  (33) 

o 

9V    =9V=9V=9V=9V  n4> 

9R      9y    9^    36  "  3<j>  l   j 


Equations  (21)  through  (34)  complete  the  formulation 

of  the  left  hand  sides  of  the  six  Lagrange  equations,  as 

generally  expressed  by  Equation  (11) .   In  completing  these 

equations  it  is  only  necessary  to  formulate  the  generalized 

forces  F„  ,  F_,  F  ,  F,,  Fn  and  Fx.       The  principle  of  vir— 
Z     R    Y    ty    8       $ 

tual  work  forms  the  basis  for  these  calculations  and  can  be 

described  mathematically  by: 
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3e 

F  (M  )  +  F  (—  )  +  F  (—  )  +  M  ( - 

XKdq    ]  YV9q  '      ZV8q  ;     X^9q 

Mr         -ir  -ir  Mr 


)  + 


My( 


3e, 

3q] 


Be, 


)  +  M„( 


Z^q-}  =  F*r 


(35) 


where  X,Y,Z  are  positive  incremental  displacements  along 

the  X,  ,Y..,Z,  axes,  e  ,e  ,e"  are  positive  incremental  anqu— 
111        xyz      ^  3 

lar  displacements  about  the  X..,Y,  ,Z,  axes. 

The  incremental  translation  or  rotation  along  or  about 
the  X,,Y   or  Z   axis,  caused  by  varying  each  generalized 
coordinate  individually,  can  easily  be  visualized  with  the 
aid  of  Figure  3.   For  example,  varying  the  generalized  co- 
ordinate Z   by  a  positive  increment  will  produce  no  trans- 
lation or  rotation  about  either  X,  or  Y, ,  and  no  rotation 
about  Z, .   Using  Equation  (35)  as  a  guide  (realizing  that 
each  X,Y  and  Z  appearing  should  be  subscripted  with  a  one, 
because  we  are  dealing  with  the  cartesian  coordinate  refer— 


!fi    !!i    !fx    ffy 

9Z   '  8Z   '  9Z   '  8Z 
o     o     o     o 


and 


ence  frame) ,  it  follows  that 

ic-=—  are  identically  zero.   From  Figure  5  it  is  clear  that 
o 


Figure  5.   Positive  Incremental  Displacement  of  the 
Generalized  Coordinate  Z  . 
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as  Z   is  incremented  a  positive  amount,  Z,  follows  direct— 
o  1 

ly  in  its  negative  defined  direction.   This  implies  that 

3Z, 

7- =  —1  and  the  generalized  force  F    =  — F    follows  from 

9Zo  o      Zl 

Equation  (35) . 

The  remaining  generalized  forces  with  respect  to  the 

cartesian  coordinate  reference  frame  (X, ,Y, ,Z,)  can  be 

derived  in  a  similar  manner.   The  complete  set  is: 

%  -  %  (36) 

Ffi    =   -M      sin    ip    +    Mv    cos    ty  (37) 

0  Xl  Yl 

F      =   M      cosQcosijj    +   M      cos6sint|j   —  M      sinQ  (38) 

9  x]_  Y^  Z^ 

Fz=-F  (39) 

o  1 

FR   =   -F^  (40) 

F      =    F      R   +    M  (41) 

Making  use  of  the  axes  transformation  [T  ,, ] ,  the  general- 
ized forces  expressed  in  terms  of  the  body  axes  are: 

FR  "  -(al2FX  +  °32FY  +  a32FZ»  (42) 

FY  =  -R(ollFx+  a21FY  +  a31FZ)+(a13MX  +  a23MY+  a33MZ>   (43) 

F*  ■  a13MX  +  "23^  +  a33MZ  (44) 

F_=  M   cos  <J>  —  M   sin  <J)  (45) 

%  "  Mx  <46> 

FZ0=  -(a13FX  +  «23FY  +  a33FZ»  <47) 
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Adding  Equations  (42)  through  (47)  to  the  respective 
right— hand  sides  of  the  earlier  mentioned  grouping  (of  the 
left— hand  terms) ,  in  the  general  expression  given  by  Equa- 
tion (11) ,  completes  the  development.   This  grouping  is 
shown  in  Appendix  A,  where  the  six  equations  of  motion  are 
presented  in  their  entirety-. 

Engine  thrust  terms  are  not  included  in  the  six  equa- 
tions of  motion.   Previous  studies  indicate  .that  use  of 
thrust  in  the  spin  is  normally  avoided  since  flameouts  are 
likely  and  serious  damage  can  result  (Ref.  4). 

D.   STEADY  SPIN  EQUATIONS 

The  full  equations  of  motion  (as  given  in  Appendix  A) 
can  be  greatly  simplified  for  a  steady  state  spin  condition. 
In  such  a  state,  the  aircraft  is  imagined  to  be  descending 
at  a  constant  velocity,  at  a  constant  radius  (R)  and  rota- 
tion rate  (y)  about  the  central  spin  axis  and  with  a  fixed 
orientation  with  respect  to  the  cartesian  coordinate  refer- 
ence frame  (X.  ,  Y.  ,  Z..  )  .   This  implies  that  all  second— deriv- 
ative  terms  and  any  first— derivative  terms  of  the  Euler 
angles  (ip,6,<J>)  can  be  set  equal  to  zero.   On  examination 
of  the  auxiliary  terms  of  Appendix  A,  it  is  easily  seen  that 
all  time  rates  of  the  direction  cosines  (a . . )  and  angular 

•      •      • 

accelerations  (p,  q,  r)  are  equal  to  zero.   After  making 
these  simplifications,  the  six  equations  of  motion  that  de- 
scribe the  steady  spin  condition  are: 

RES(l)  =  0  =  a13Fx  +  a23Fy  +  «33FZ  -  mg  (48) 
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RES(2)     =    0    =    ot12Fx    +    ct22Fy    +    a32Fz    -  my2R  (49) 

RES (3)    =    0   =    R(a11Fx+a21Fy    +    a31FZ)     + 

(a13MX    +a23MY    +    a33MZ}  (50) 

RES  (4)    =    0    =ct13Mx   +    a23MY    +    a33Mz)  (51) 

RES  (5)     =    0    =   -Y(Ivpcos0"   +    I    qsin<j>sin6    + 

I    (rcoscj)sinG)    +   M   costp   —  M   sin<J)  (52) 

Li  X  £j 

— Ylvl7  (pcosc{)sin0  +  rcosG) 
xz 

RES  (6)  =  0  =  Mx-qr(Iz-Iy)-IxzPYsin4>cos9  (53) 

The  aerodynamic  coefficients  are  written  in  a  form  con- 
sistent with  the  data  obtained  from  Ref.  5.  The  aerodynam- 
ic forces  and  moments  are  written: 

Fx  =  |PV2S(CX  +  C    •  6e  +  Cxq  •  q)  (54) 

Fy  =  ipV2S(Cy+C    -5  +C    -«r+C   -p  +  Cy  -r)  (55) 

o        o        p       r 

Fz  =  ^PV2S(CZ  +  C    -fi   +  Cz  -q)  (56) 

e        2 

MX  =  2lpv2sb(C£+  C£   *6r+  C£*P  +  Cl    'r)  (57) 

■  <$r        P       r 

^  =  ^-pV2Sc(C  +  C    «6   +  C   -q)  (58) 

i         2  m   m .    e    m   ^ 

6e        q 

M_  =  ^pV2Sb(C  +C    -6  +C    -6  +C   • p+C   -r)  (59) 

Z    2       n   n~    a   nr    r   n  e      n 

6        6        p      r 
a       r      * 

If  the  relations  given  by  Equations  (54)  through  (59) 
were  substituted  into  the  steady  spin  equations,  as  given 
by  Equations  (48)  through  (53),  the  six  equations  would 
then  describe  the  resulting  aircraft  state  in  terms  of  con- 
ventional aerodynamic  forces  and  moments. 
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The  residual  notation  (RES(i))  introduced  in  the  steady 
spin  equations,  as  given  by  Equations  (48)  through  (53),  is 
used  as  a  measure  of  how  closely  each  associated  equation 
is  being  satisfied.  Ideally,  each  residual  should  be  iden- 
tically equal  to  zero  but  this  is  not  possible  due  to  trun- 
cation and  round— of f  error  within  the  computer  at  various 
stages  of  the  program. 

A  criteria  function  is  defined  that  will  indicate  how 
close  the  complete  set  of  six  equations  is  to  a  steady  spin 
condition.   The  criteria  function  chosen  for  this  investi- 
gation was  the  sum  of  the  squares  of  the  individual  resid- 
uals, which  is  written: 

F  =  RES(1)2+  RES(2)2+  RES(3)2+  RES(4)2+  RES(5)2  +  RES(6)2 

(60) 
This  function  will  be  the  object  of  the  minimization  effort 
in  the  optimization  procedure  described  in  the  next  section 

The  airplane  has  six  degrees  of  freedom  with  second 
order  differential  equations;  therefore,  12  variables  are 
required  to  specify  its  state.   For  convenience  the  12 

•        •        •  •  •  • 

variables  chosen  are:  \p  ,\p  ,Q  ,Q  ,<$>,$  ,Z    ,  Z  ,R,R,  y  and  y.      Mak- 
ing use  of  the  fact  that  each  equality  constraint  imposed 
upon  the  system  will  remove  one  variable,  the  total  number 
of  required  variables  can  therefore  be  reduced.   Some  of 
these  have  been  discussed  in  describing  the  steady  spin; 
specifically,  that  R=if;  =  0  =  <j>  =  O.   They  were  also  used 
to  reduce  the  general  equations  given  in  Appendix  A  to  the 
steady  spin  equations  given  by  Equations  (48)  through  (53) . 
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By  specifying  that  the  search  be  made  at  y    =  0  and  at  a 
specified  altitude  (Z  ) ,  the  total  number  of  equality  con- 
straints is  six.   This  leaves  six  variables  to  search  over 
if  each  control  surface  is  held  in  a  fixed  position.   The 
chosen  set  of  variables  to  be  searched  over  in  this  inves— 
tigation  are:  y,    R,  Z  ,  \\)J   Q    and  <j) .   Specifying  these  six 
initial  conditions,  the  altitude  parameter  (p,  as  Z 
doesn't  appear  explicitly  in  the  steady  spin  equations), 
the  physical  parameters  of  the  aircraft  (mass,  chord,  span, 
etc.),  and  having  a  set  of  tabulated  aerodynamic  coeffi- 
cients, the  criteria  function  can  be  evaluated.   The  only 
thing  now  required  is  a  systematic  method  of  varying  the 
six  variables  (y ,  R,  Z  ,  \p ,  0,  4> )  to  obtain  the  minimum 
value  of  the  criteria  function  close  to  zero. 

The  solution  of  the  steady  spin  condition  involves  the 
simultaneous  solution  of  six  highly  coupled  nonlinear  dif- 
ferential equations  involving  aerodynamic  data.   There 
does  not  exist  a  closed  form  analytic  solution  to  this  set 
of  equations  and  therefore  a  numerical  method  must  be 
employed. 
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III.   REVISED  FORMULATION  OF  THE  STEADY  SPIN  EQUATIONS 

The  steady  state  equations  were  FORTRAN  coded  for  use 
in  program  SPIN  and  initial  efforts  at  predicting  the  equi- 
librium spin  condition  were  unsuccessful.   During  these 
efforts  it  was  realized  that  the  rate  of  convergence  was 
extremely  slow  and  that  the  solutions  were  not  close  to  the 
states  already  known  to  exist,  as  shown  by  Adams  (Ref.  7). 

After  some  analysis  it  was  theorized  that  the  deletion 
of  the  acceleration  terms  from  the  derived  equations  of 
motion  was  the  cause.   The  equations  should  have  been  un- 
coupled and  solved  explicitly  for  these  acceleration  terms. 
Then  the  residual  terms,  representing  the  acceleration 
values,  would  be  a  true  measure  of  just  how  close  the  air- 
craft is  to  the  equilibrium  spin  condition. 

Attempts  at  uncoupling  the  full  equations  of  motion,  as 
given  in  Appendix  A,  showed  that  it  would  be,  at  most,  a 
very  time  consuming  and  nearly  impossible  task  as  the  \p ,    9 
and  4>  equations  are  so  highly  coupled.   Instead,  a  search 
for  an  alternate  set  of  moment  equations  was  made  with  the 
decision  to  use  the  equations  given  by  Etkin  (Ref.  3).   The 
three  equations  express  the  moments  about  the  body  axes  and 
are : 


MX  =  XXP-  XXZr  +  W  V  -  ^W  (6D 
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My  =  Iyq  +  rp(Ix-  Iz)  +  IXZ(P2~  q2)  (62) 

MZ  =-ZXzP  +  V  +    Pq(IY~  V  +  ^r  (63» 

The  first,  three  equations  of  the  original  formulation 
express  Newton's  second  law  of  motion  along  the  radial, 
tangential  and  vertical  directions  through  the  aircraft's 
center  of  gravity.   These  equations  can  be  expressed, 
without  introducing  the  transformations  (a. .),    about  the 
cartesian  coordinate  reference  frame  as: 

Fv   =  my  (64) 

Xl 

•  2 
F    +  my  R  =  mR  (65) 

1 

F    =  mZ  (66) 

Zl 

These  three  force  equations  when  combined  with  the  three 
moment  equations,  as  given  by  Equations  (61)  through  (63), 
describe  the  motion  of  and  rotation  about  the  center  of 
gravity  of  the  aircraft  and  therefore  completely  describe 
the  motion  of  the  six— degree— of— freedom  vehicle. 

By  explicitly  solving  Equation  (64)  for  the  accelera- 
tion term  y,  Equation  (65)  for  R,  Equation  (66)  for  Z  and 
Equation  (61)  through  (63)  for  p,  q  and  r  the  equations 
are  in  the  form  desired  for  the  optimization  search. 
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IV.   OPTIMIZATION  PROCEDURE 

A.   CONCEPT  OF  OPTIMIZATION 

Using  state  variable  notation,  the  equilibrium  solution 
and  optimization  scheme  f or  this  particular  formulation  can 
be  described  mathematically  as  follows.   If  we  let 


/!\ 


X 


-<    - 


R 
Z 
P 

q 


> 


,  the  acceleration  terms 


and 


I  *  \ 


-< 


R 
Z 


,    the  independent  variables , 


\   * 

therefore    X    =  ]  F (Y) [   .   Then  the  optimization  scheme 
finds  ]  Y  (  such  that  1  X  [  =0,  for  equilibrium.   The  op- 
timization scheme  employed  in  this  investigation  uses  the 
criteria  function  which  is  the  sum  of  the  squares  of  the 

acceleration  terms  (modified  equation  residuals) ,  given  by 

6 

x.   =  y  +  R  +  Zz+  p  +  q  +  r 
i=l    1 
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and  finds  )    Y  V  such  that  J (Y)  is  an  absolute  minimum 
which  yields    X    =0. 


B.   OPTIMIZATION  ALGORITHM 

The  optimization  method  used  in  this  investigation 
finds  the  minimum  value  of  the  criteria  function,  which  is 
a  real— valued  function  of  several  variables.   The  criteria 
function  is  used  as  an  indicator  of  how  close  the  aircraft 
is  to  being  in  an  equilibrium  spin  condition.   The  minimi- 
zation process  takes  into  account  any  constraints  imposed 
upon  the  system,  which,  in  this  case,  were  imposed  limits 
of  angle  of  attack  and  sideslip  from  the  available  aero- 
dynamic coefficients. 

The  optimization  procedure  starts  with  a  given  set  of 
independent  control  parameters  and  varies  the  current  values 
of  these  parameters,  within  the  imposed  boundaries,  to  im- 
prove the  value  of  the  criteria  function.   The  process  of 
parameter  adjustment  continues  until  the  minimum  of  the 
criteria  function  is  reached.   Depending  on  the  hyper— sur- 
face defined  by  the  independent  control  parameters,  the 
minimum  found  may  be  just  a  local  minimum  which  would  not 
satisfy  the  desired  steady— spin  criteria  functional  value 
(equal  to  zero) .   The  absolute  minimum  of  the  hyper— surf ace 
will  yield  a  criteria  function  value  equal  to  zero  and  is 
the  object  of  the  search.   The  initial  conditions  given 
for  the  independent  control  parameters  will  greatly  influ- 
ence the  success  of  finding  the  absolute  minimum.   There- 
fore, there  is  a  certain  trial  and  error  process  in  the 
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search  for  the  absolute  minimum  without  any  guarantee  that 
one  exists.   This  then  is  the  main  disadvantage  of  the  op- 
timization procedure  because  in  spite  of  the  energy  and 
perseverance  of  the  researcher  he  cannot,  at  the  end,  say 
that  a  steady  spin  condition  does  not  exist  for  that  air- 
craft if  in  fact  he  hasn't  been  able  to  locate  an  absolute 
minimum.  . This  is  a  very  common  disadvantage  shared  by  all 
optimization  procedures  available  at  this  time. 

The  optimization  algorithm  used  in  the  spin  program  is 
called  EXTREM.   A  complete  listing  of  this  subroutine  is 
given  in  Appendix  B  along  with  a  functional  flow  chart.  The 
user  specifies  a  set  of  independent  variables  in  name  and 
initial  value  and  a  corresponding  number  of  scaled  incre- 
ments defining  the  initial  step  sizes  along  each  of  the  in- 
dependent variables.   The  algorithm  first  checks  that  the 
given  initial  conditions  do  not  violate  any  constraints 
imposed  on  the  system  and  then  proceeds  to  determine  the 
main  line  of  search  in  the  variable  space.   The  initial  con- 
ditions form  a  point  in  this  variable  space  and  the  speci- 
fied initial  step  sizes  give  a  second.   A  vector  joining 
these  two  points  determines  the  initial  main  line  of  search 
for  the  first  stage,  as  shown  in  Figure  6  for  a  two— variable 
problem.   Along  this  main  line  the  approximate  minimum  point 
is  found  by  extrapolating  a  parabola  through  three  points 
determined  by  X. -  DX,  X,  and  X,+  DX  and  then  interpolating 
for  the  minimum  (X.  ,).   A  Gram— Schmidt  orthogonalization- 
process  determines  the  secondary  direction  of  the  first 
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staqe  which  is  orthogonal  to  the  main  direction  at  X .  n . 
The  approximate  minimum  point  X .  ~  is  found  along  this  line 
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Figure  6.   Main  Line  of  Search  in  a  Two  Variable  Problem 

For  a  third  variable  (not  shown  in  Figure  6)  a  line  of 
search  through  X .  „  is  determined  by  the  Gram— Schmidt  pro- 
cess that  is  orthogonal  to  both  the  main  and  secondary 
directions.   Each  new  direction  is  orthogonal  to  all  the 
preceding  directions  of  search  and  the  procedure  for  deter- 
mining these  lines  of  search  and  approximate  minimum  points 
is  continued  until  the  minimum  along  the  K    line  (K  being 
the  number  of  independent  variables)  is  found.   At  this 
point  the  first  stage  is  complete  and  the  main  direction 
for  the  second  stage  is  determined  by  a  line  joining  the 

point  X,  and  X.  ,,  ,    where  X.  ,,  is  the  approximate  minimum 
r      1      i+k         l+k         ^ 

along  the  K    line  of  the  first  stage.   The  second  stage 
proceeds  just  as  described  for  the  first,  and  the  overall 
stage  procedure  continues  until  a  local  or  absolute  rela- 
tive minimum  is  found.   Between  stages  the  program  checks 
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the  stopping  conditions  on  the  criteria  function  value, 
the  arguments,  and  the  maximum  number  of  stages  allowed  to 
reach  the  minimum,  and  will  stop  if  any  of  these  conditions 
is  met.   The  step  size  along  a  line  of  search  is  variable 
in  the  program,  which  allows  the  search  steps  to  decrease 
near  the  minimum  or  within  tight  curves  and  to  increase 
again  once  these  curves  are  past. 

A  more  detailed  description  of  the  optimization  method, 
the  algorithm  and  several  test  cases  can  be  found  in  Ref.  8 
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V.   DISCUSSION  OF  RESULTS 

A.   VERIFICATION  EFFORT 

In  an  attempt  to  verify  that  the  spin  equations  of 
motion  describe  the  aircraft  motion  and  also  to  check  out 
the  entire  computer  program,  several  runs  were  made  for  an 
equilibrium  flight  condition  of  straight  and  level  flight. 
The  choice  of  this  flight  mode  was  made  so  that  the  results 
could  be  verified  by  hand  while  searching  over  a  reduced 
number  of  independent  variables  in  satisfying  the  same  cri- 
teria function  as  discussed  earlier. 

As  mentioned  earlier,  thrust  terms  were  not  included  in 
the  general  equations  of  motion,  as  previous  studies  indi- 
cated that  its  use  is  normally  avoided,  since  flameouts  are 
likely  and  serious . engine  damage  can  result.   However,  in 
order  to  simulate  straight  and  level  flight  thrust  was  in- 
cluded with  the  very  restrictive  assumption  that  it  was 
orientated  in  the  negative  X  direction  and  passed  through 
the  center  of  gravity. 

Straight  and  level  flight  can  be  simulated  in  program 

o         • 
SPIN  by  specifying  a  value  of  rp    equal  to  90  ,  y    and  Z  both 

set  to  zero.   By  assuming  that  a  wings— level,  straight  line- 

of— flight  orientation  can  be  obtained  with  neutral  rudder 

and  aileron  control,  the  value  of  6   and  6   are  set  equal 

a      r  ^ 

to  zero  as  well  as  the  bank  angle  <p  .   The  value  of  the 
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radius  vector  (R)  is  arbitrary.   With  this  accomplished  for 
each  run  a  value  of  the  pitch  angle  (G)  (in  this  special 
case  0  is  the  same  as  the  angle  of  attack)  is  specified. 
Three  independent  variables  (thrust  required,  elevator 
angle  and  airspeed)  and  their  initial  step  sizes  are  speci- 
fied in  the  main  program  and  the  routine  yields  the  value 
of  these  three  variables  that  would  be  required  for  straight 
and  level  equilibrium  flight. 

The  results  of  these  runs  were  verified  by  using  an 
electronic  calculator;  the  equations  of  motion  were  satis- 
fied.  The  numerical  results  of  thrust,  velocity  and  eleva- 
tor deflection  angle  were  reasonable  for  that  particular 
flight  configuration. 

B.   PROGRAM  SPIN  RESULTS 

Several  SPIN  runs  were  conducted  for  configuration  B  in 
order  to  compare  the  results  with  those  given  by  Adams 
(Ref.  9).   The  mass,  dimensional  characteristics  and  aero- 
dynamic coefficients  for  this  aircraft  configuration  are 
listed  in  Appendix  C. 

The  results  of  the  predicted  spin  characteristics  for 
configuration  B  are  shown  in  Table  I  where  they  are  com- 
pared with  the  results  obtained  by  Adams.   The  predicted 
spin  characteristics  developed  by  program  SPIN,  using  the 
optimization  algorithm  EXTREM,  are  in  good  agreement  with 
those  predicted  by  Adams.   The  slight  difference  in  the 
first,  second  or  third  decimal  digit  in  some  cases  is  prob- 
ably  a  result  of  the  particular  computer  being  used. 
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Adams  used  a  CDC  machine  with  a  bigger  bit  word  than  the 
IBM— 360  and  therefore  more  significant  figures  could  be 
handled  in  the  numerous  calculations  performed.   As  a  re- 
sult of  this  the  minimum  value  of  the  criteria  function 

—24 
obtained  with  program  SPIN  was  10    .   Adams  has  been  able 

to  drive  his  criteria  function  residuals  (also  a  sum  of 

the  squares  but  of  a  different  set  of  terms)  to  a  value  of 

io-36. 

The  third  grouping  of  results  in  Table  I  show  very 
little  similarity,  despite  the  same  control  settings  being 
used.   The  first  steady  spin  predicted  is  an  entirely  dif- 
ferent equilibrium  condition  as  can  be  seen  by  comparing 
the  values  of  the  radius  (R)  and  angular  rate  (y) .   The 
second  spin  predicted  would  at  first  appear  to  be  the  same 
as  that  given  by  Adams.   The  difference  in  values  of  the 
angle  of  attack  (a) ,  inertial  velocity  (V)  and  radius  (R) 
make  it  an  entirely  different  spin  condition  and  not  just 
a  poor  representation  of  the  same  spin.   This  statement  is 
justified  by  the  closer  correlation  of  results,  in  the 

other  groupings,  for  the  same  value  of  the  criteria  func— 

—24 
tion  (10    )  . 

The  computer  prediction  runs  for  configuration  A  are 

shown  in  Table  II.   There  were  no  other  results  to  compare 

with  this  configuration.   The  spin  prediction  scheme  was 

applied  to  this  particular  aircraft  to  assure  that  the 

method  would  at  least  work  for  another  aircraft.   The  choice 
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for  this  second  configuration,  over  other  possibilities, 
was  made  easy  because  the  aerodynamic  derivatives  were  al- 
ready in  tabulated  form  as  given  in  Ref.  7. 

The  effects  of  density  variation  on  the  equilibrium 
spin  conditions  is  shown  in  Table  III.   Grantham  and  Graf- 
ton (Ref.  8)  concluded  from  their  studies  that  a  decrease 
in  density  gave  a  faster  rotation  rate  (y) ,    higher  veloc- 
ities and  little  change  in  angle  of  attack  and  sideslip. 
The  parameters  predicted  by  program  SPIN  verify  this  be- 
havior with  the  added  effect  that  the  radius  of  the  spin 
tends  to  decrease,  i.e.,  the  spin  becomes  tighter  with  a 
decrease  in  density.   The  change  in  the  orientation  angle 
0  would  indicate  that  the  spin  would  be  somewhat  flatter 
with  a  decrease  in  density. 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

As  a  result  of  developing  and  presenting  the  analytic 
spin  prediction  technique,  the  following  conclusions  and 
observations  are  made. 

1.  The  original  developed  equations  of  motion,  as 
presented  in  Appendix  A,  are  too  highly  coupled 

to  be  an  effective  tool  in  this  type  of  prediction 
scheme.   In  this  light,  those  equations  can  be  con- 
sidered, at  most,  as  an  academic  exercise  leading 
to  a  rather  unique  formulation. 

2.  The  comparison  of  results  demonstrated  the  accu- 
racy of  the  method  and  the  effects  of  density 
variation  study  showed  one  possible  use  of  the  pro- 
gram.  The  capability  of  this  method  to  predict 
steady  spin  conditions  is  limited  only  by  the  accu- 
racy with  which  the  mass,  physical  characteristics 
and  aerodynamic  coefficients  are  determined. 

3.  Due  to  the  direct  modeling  in  a  cylindrical  coordi- 
nate reference  frame  this  method  has  a  certain 
advantage  over  that  demonstrated  by  Adams  (Ref.  7). 
The  number  of  equality  constraints  on  the  solution 
required  by  this  prediction  method  is  far  fewer, 
which  allows  for  a  much  simpler  computer  formula- 
tion and  equally  accurate  results.   However,  this 
method  uses  more  computer  time  than  Adams '  in 
converging  to  the  equilibrium  spin  condition. 

4.  As  with  any  optimization  technique  the  importance 
of  equally  weighing  the  residuals  of  the  equations 
is  of  utmost  importance.   The  value  of  the  residuals 
should  be  as  close  as  possible  to  each  other,  in 
order  of  magnitude,  for  proper  convergence  of  the 
optimization  algorithm.   Unequal  weighting  of  any 

of  these  residuals  cause  these  terms  to  dominate 
and  influence  the  search  in  a  direction  that  doesn't 
actually  contain  the  desired  solution  minimum.  Many 
hours  of  research  time  were  lost  by  not  recognizing 
this  important  fact. 

In  light  of  the  demonstrated  usefulness  of  the  analytic 

prediction  scheme  the  following  follow— on  work  could  prove 
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valuable  in  better  understanding  the  spin  problem.   It  is 
recommended  that: 

1.  The  use  of  rotary— balance  aerodynamic  data  be  in- 
vestigated as  a  possible  way  of  more  closely  relat- 
ing the  aerodynamic  coefficients  of  the  tunnel 
model  to  that  experienced  by  the  real  airplane. 

The  use  of  this  type  of  data  will  allow  at— desk 
calculation  of  the  steady  spin  conditions.   This 
particular  spin  estimation  method,  if  used  to  pre- 
dict the  initial  conditions  for  program  SPIN, 
would  reduce  any  sensitivity  to  the  guessed  initial 
conditions  and  also  reduce  the  required  computer 
convergence  time.   This  particular  method  of  spin 
estimation  is  given  in  Ref.  9. 

2.  The  equations  of  motion  could  be  linearized  in 
order  to  obtain  some  information  about  the  stability 
of  the  predicted  equilibrium  spin  condition. 

3.  This  prediction  scheme  could  be  coupled  with  exist- 
ing fixed— base  spin  simulators  to  provide  the 
starting  initial  conditions  for  the  spin  simulation. 
With  this  the  recovery  characteristics  could  be 
investigated,  recovery  techniques  perfected  and  the 
necessary  pilot  interactions  evaluated  and  studied. 
This  one  tool  alone  could  provide  valuable  insight 
into  the  problem  areas  and  shows  much  promise  in 
future  work. 
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APPENDIX  A 

FULL  EQUATIONS  OF  MOTION 

The  following  six  equations  completely  describe  the 
rigid  body  motions  of  an  aircraft,  using  body  axes,  as 
formulated  in  a  cylindrical  coordinate  reference  frame. 

Z  EQUATION 

mZo  =  Q13FX  +  a23FY  +  a33FZ  +  mg 
R  EQUATION 

mR  -  mY2R  =  -(a^  +  a22Fy  +  a^F^ 

y  EQUATION. 

mR(Ry  +  2Ry)  +  IX(P«13  +  Pa13)  +  IY  (q(*23  +  qa23) 

•  •         • 

+  I   {— p(0cos<j)sin9  +  <J)sinc{)COs6 )  + 
p  cost|)cos0  —  9rcos9  —  rsin6}  = 

R(allFX  +  a21FY  +  a31FZ}  +  a13MX  +  a23MY  +  a33MZ 
ip    EQUATION 

Ix(Pa13  +  pa13)  +  Iy(qa23  +  qa^)  + 

Iz(ra33  +  ra33)  +  I   {-p  (9coscJ>sin9  + 

(J>sin(})COs9)  +  p  cosijicosG  —  9r  cosG  —  r  sin9}   = 
a13Mx  +  a23My  +  a33Mz 
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6    EQUATION 


I     (q    coscj)   —  q<{>sin<{>)    —   I     (rsincf)    +    rcfjcosc})) 

Y  LA 

+    (i'+Y)  {Ivpcos0    +    I    qsin<J)sinQ    +    I    rcos<J)sinG} 
-Ixz  (<}>cos(J)p   +    sin^p)    +    I       (pcos<J)sin9  (i+y) 
+    rcosG  (ip+y)  )    =   M      coscj)   —  M      sintf) 


4)    EQUATION 

IXP    +    5r(Iz    ~    IY)    _   JXZ(^    +    P^00^    +    sin<f>cos6  (i+y)  )  ) 

=    M 
X 

where, 

a,  ,     =  COSij;COS9 

a,  ,    =  — i|>sirnjjcos0   —    9cosiJjsin9 

a,  2   =  sinijjcosG 

•  •  • 

a,  ~    =  i|;cosipcos9   —    0sintjjsin9 

a,  -.    =  — sin9 

a, 3   =  — 9cos9 

a~,    =  cos^sinGsincf)   —   sini/jcoscf) 

•  •  • 

a?n    =   — ^sinipsinGsincJ)   —    9cos^cos9sincf) 

•  •  • 

— 4>cos^sin9cosi[i  —   ipcosi|;cos(f>    +    4>sirnpsin<J) 
a99    =    sini|;sin0sin'(J)+   cosijjcoscj) 

•  •  • 

a99    =    ^cos^sin9sin(|)—   9sinipcos9sin({) 

•  •  • 

—   4>siniJ;sin9cos(f)   —   ^sin^cos^   —   ^cos^sincf) 
a--.    =   cos9sin(J) 
a9^   =    cf>cos9cos<f>   —   9sin9sin<j) 
a.,,    =    cos^sinGcoscj)    +    sinijjsin(f> 

•  •  • 

a-.,    =  — ipsiruJjsinQcoscf)    +    GcosijjcosGcoscJ) 

•  •  • 

— <|)COSi|jsin9sin(J)    +    i^cos^sin^    +    ^3inii;cos<{> 
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a 


a 


a 


a 


32 

32 

33 

33 

P 

P 

q 
q 


r 

r 


sin^sinGcoscf)   —      cosijjsincf) 

^cosipsinGcos^    +    OsimpcosGcoscJ)   —   4>sin(f)sin6sin< 

+    ^sinipsincf)   —   ^cosipcoscf) 

cos9cos<j) 

—  GsinGcoscf)   —   <|)sin(})cos0 
4>   —   sinG  {ty+y) 

<f>   —   ipQcosG    —   ijjsinG   —   y9cos6    —  ysinQ 

•  •         • 

GcoscJ)    +      cosGsin<|)  (ip+y) 

—  6<f)sincf>    +    0costj)   —   ^GsinGsincf)   —      ^cosQcosc}) 
+    ^cosGsint})   —  yGsintfisinG    +      y^cosQc03^ 

+    YCOsQsinc}) 

•         •  • 

cosQ    cos<j>  (ty+y)    —   GsincJ) 

—  iJjQsinGcoscf)   —   i^cosOsirn})    +    ^cosGcoscf) 

—  Oncost})   —   Osintj)   —      YQsi11^003^ 

—  Y^cosOsin^    +      y005©003^ 
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APPENDIX  B 


PROGRAM  SPIN 


A.  INTRODUCTION 

The  intent  of  this  appendix  is  to  serve  as  a  users 
manual  or  guide  to  the  use  of  Program  SPIN.   An  alphabeti- 
cal listing  of  the  computer  variables  and  their  meanings 
is  provided.   A  brief  description  of  the  main  purpose  of 
each  subprogram  as  well  as  an  extensive  explanation  on  pre- 
paring the  data  deck  is  included,   A  listing  of  the  program 
with  numerous  comments  is  included  along  with  several  func- 
tional flow  charts  to  help  in  understanding  the  structure 
and  relationships  of  the  subroutines.   A  sample  output  is 
given  and  the  storage  and  time  requirements  of  the  program 
are  discussed .  . 

B.  LIST  OF  COMPUTER  VARIABLES 

The  following  is  an  alphabetical  listing  of  the  compu- 
ter variables  and  their  meaning. 
All— A33      Direction  cosines 

A(  )         Working  matrix  for  Subroutine  EXTREM 
ALPHAD       Aircraft  angle  of  attack,  degrees 
ALPHAR      Aircraft  angle  of  attack,  radians 
B  Wing  span,  ft 

BETAD        Angle  of  sideslip,  degrees 
BETAR        Angle  of  sideslip,  radians 
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CBAR  Wing  chord,  ft 

CL  Total  roll  moment  coefficient 

CM  Total  pitch  moment  coefficient 

CN  Total  yaw  moment  coefficient 

COEFA (  )     Aerodynamic  coefficients  which  are  a  function 
of  angle  of  attack  only 

COEFAB (  )    Aerodynamic  coefficients  which  are  a  function 
of  angle  of  attack  and  angle  of  sideslip 

CPHI  Cosine  of  PHI 

CPSI  Cosine  of  PSI 

CR(  )  Working  vector  of  interpolated  coefficients 

CTIIETA  Cosine  of  THETA 

CX  Total  X  force  coefficient 

CY  Total  Y  force  coefficient 

CZ  Total  Z  force  coefficient 

D  1.0-(ZXI)**2/XI*ZI 

DDELA        Incremental  step  size  of  aileron  control,  de- 
grees 

DDELE        Incremental  step  size  of  elevator  control, 
degrees 

DDELR        Incremental  step  size  of  rudder  control,  degrees 

DELTAA      Aileron  control  deflection,  degrees  (Positive 
with  right  aileron  trailing  edge  down) 

DELTAE       Elevator  control  deflection,  degrees (Positive 
with  trailing  edge  down) 

DELTAR       Rudder  control  deflection,  degrees  (Positive 

with  trailing  edge  left  when  viewed  from  above) 

DFMAX        Stopping  condition  of  the  functional  variation 
in  Subroutine  EXTREM 

DGAMD        Incremental  step  size  of  GAMDOT,  rad/sec 
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DPHI 
DPS  I 

DR 

DRDOT 

DTHETA 

DTHRUS 

DXMAX 

DX(  ) 

DZDOT 

F 

FMAX 

FOPT 
G 

GAMDOT 

HRHOS 
IW 


K 

LMAX 

NPRINT 

P 

PBAR 


Incremental  step  size  of  PHI,  degrees 

Incremental  step  size  of  PSI,  degrees 

Incremental  step  size  of  R,  ft 

Incremental  step  size  of  RDOT,  ft/sec 

Incremental  step  size  of  THETA,  degrees 

Incremental  step  size  of  THRUST,  lbs 

Control  parameter  for  Subroutine  EXTREM 

Array  of  step  sizes  of  the  independent  variables 

Incremental  step  size  of  ZDOT,  ft/sec 

Function  to  be  minimized 

Desired  minimum  value  of  the  optimization 
function 

Optimum  value  of  the  function  being  minimized 

Gravitational  acceleration,  ft/sec2 

Rate  of  change  of  the  cylindrical  orientation 
variable  y,    rads/sec 

0.5*RHO*S 

Control  parameter  of  Subroutine  EXTREM  for 

writing  instructions 

+1   All  output  suppressed 

+2   Final  output  only 

+3   Output  at  the  end  of  each  stage 

The  number  of  independent  variables  being 
searched  over 

Stopping  condition  on  the  number  of  stages  of 
Subroutine  EXTREM  (A  negative  sign  indicates 
that  a  minimum  is  sought) 

Print  indicator,  set  in  MAIN  and  passed  to 
Subroutine  DATAIN  to  echo  check  the  aerodynamic 
coefficients  (  =0,  No  print;  =1,  Prints) 

Angular  velocity  about  the  X  body  axis,  rad/sec 

P*B/2V 
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PHI 
PHIR 
PSI 
PSIR 

Q 

QBAR 

R 

RBAR 

RDOT 

RES(  ) 

RHO 

RZ 

S 

SPHI 

SPSI 

SSAM 

STHETA 

SUMRES 

THETA 

THETAR 

V 

VSQ 

VXB 

VYB 


Euler  orientation  angle  <\> ,  degrees 

Euler  orientation  angle  §,  radians 

Euler  orientation  angle  \p ,  degrees 

Euler  orientation  angle  ip ,  radians 

Angular  velocity  about  the  Y  body  axis,  rad/sec 

—  i 

Q*CBAR/2*V 

Cylindrical  radius,  ft 

R*B/2*V 

Rate  of  change  of  the  cylindrical  radius, 
ft/sec . 

Residual  of  any  of  the  spin  equations 

2    4 
Density  of  air,  lbs— sec  /ft 

Angular  velocity  about  the  Z  body  axis,  rad/sec 

2 
Wing  area,  ft 

Sine  of  PHI 

Sine  of  PSI 

Aircraft  mass,  slugs 

Sine  of  THETA 

Sum  of  the  squares  of  the  individual  equation 
residuals 

Euler  orientation  angle  6 ,  degrees 

Euler  orientation  angle  6 ,  radians 

Velocity  of  the  vehicle  center  of  mass,  ft/sec 

Velocity  squared 

Vehicle  velocity  component  along  the  X  body 
axis,  ft/sec 

Vehicle  velocity  component  along  the  Y  body 
axis,  ft/sec 
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VYVB 
VZB 

VZXB 
X 

XFORCE 
XI 

XMOM 

YFORCE 
YI 

YMOM 

ZDOT 

ZFORCE 
ZI 

ZMOM 

ZXI 


VYB/V 

Vehicle  velocity  component  along  the  Z  body 
axis,  ft/sec 

VZB/VXB 

Working  vector  that  contains  the  independent 
variables 

Aerodynamic  force  along  the  X  body  axis,  lbs 

Moment  of  inertia  about  the  Z  body  axis, 
slugs— ft2 

Aerodynamic  moment  about  the  X  body  axis, 
ft-lbs 

Aerodynamic  force  along  the  Y  body  axis,  lbs 

Moment  of  inertia  about  the  Y  body  axis, 
slugs— ft2 

Aerodynamic  moment  about  the  Y  body  axis, 
ft-lbs 

Vehicle  velocity  in  the  cylindrical  axis 
direction,  ft/sec 

Aerodynamic  force  along  the  Z  body  axis,  lbs 

Moment  of  inertia  about  the  Z  body  axis, 
slugs— ft2 

Aerodynamic  moment  about  the  Z  body  axis, 
ft-lbs 

Cross  product  of  inertia  about  the  X  and  Z 
body  axes,  slugs— ft2 


C.   PROGRAM  DESCRIPTION 

Program  SPIN  consists  of  a  main  controlling  program  and 
six  subprograms  as  illustrated  schematically  in  Figure  Bl. 

The  MAIN  program's  most  important  function  is  to  set 
the  independent  search  variables  and  their  initial  step  size 
Otherwise  the  MAIN  program  is  basically  an    input/output 
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routine.   The  initial  conditions  and  optimization  control 
parameters  are  read  into  common  memory.   The  initial  con- 
ditions are  printed  out  as  are  various  final  conditions 
of  the  airplane. 

Subroutine  DATAIN  is  called  by  the  MAIN  program  to  read 
in  the  aerodynamic  coefficients.   Coefficients  that  are  a 
function  of  both  alpha  and  beta  are  read  into  the  COEFAB 
array  and  those  that  are  a  function  of  alpha  only  are  read 
into  the  COEFA  array.   DATAIN  will  write  out  all  of  the 
coefficients  if  the  variable  NPRINT  is  set  equal  to  one. 

Subroutine  EXTREM  is  the  optimization  program  which 
determines  the  minimum  of  the  criteria  function  without 
having  to  calculate  any  derivatives.   The  single  most  im- 
portant function  of  this  subroutine  is  in  determining  the 
optimum  direction  of  search  since  the  criteria  function 
being  minimized  resides,  and  is  evaluated,  in  Subroutine 
FN.   The  returned  values  of  the  criteria  function  are  used 
to  interpolate  a  parabola  along  a  line  of  search,  to 
determine  the  extremum  on  this  parabola  and  to  determine 
new  directions  of  search. 

Subroutine  FN  evaluates  the  criteria  function  value 
through  calls  to  Subroutine  CONST,  Subroutine  CLOPCUP  and 
Subroutine  SSEQNS.   In  addition,  boundaries  on  the  search 
are  -imposed  due  to  the  limited  range  of  angle  of  attack 
and  sideslip  angle  aerodynamic  coefficients  available. 

Subroutine  CONST  calculates  various  equation  parameter 
values  needed  in  the  equations  of  motion.   Values  of 
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angle  of  attack  and  sideslip  are  determined  so  that  the 
limits  of  the  available  coefficients  aren't  exceeded. 

Subroutine  CLOKUP  takes  the  passed  values  of  angle  of 
attack  and  sideslip  and  performs  a  linear  interpolation 
on  all  the  aerodynamic  coefficients. 

Subroutine  SSEQNS  contains  the  six"  equations  of  motion 
of  the  airplane  equated  to  the  residual  (RES (I),  1=1,6) 
vector.   The  residuals  are  evaluated  and  passed,  through 
common  memory,  to  Subroutine  FN  so  that  the  criteria  func- 
tion can  be  evaluated. 

Functional  flow  charts  of  the  MAIN  program  and  each 
subroutine  are  shown  in  Figure  B2  through  B8 . 


55 


MAIN 

PROGRAM 
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SUBROUTINE 
DATA IN 
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PRINT 
OUT_- 
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COEFFICIENT 
PRINT  OUT 

SUBROUTINE 
EXTREM 

' 

: 

SUBROUTINE 
FN 

] 
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1 

SUBROUTINE 
CONST 

SUBROUTINE 
CLOKUP 

SUBROUTINE 
SSEQNS 

Figure  Bl.   Schematic  of  Program  SPIN 
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START 


INITIALIZE  SEARCH  VARIABLE  STEP  INCREMENTS 


READ  INITIAL  CONDITIONS  AND 
EXTREM  CONTROL  PARAMETERS 


SET  COEFFICIENT  PRINT  INDICATOR 
SET  THE  NUMBER  AND  SPECIFY  THE 
SEARCH  VARIABLES 


CALL  SUBROUTINE  DATAIN 


PRINT  THE  INITIAL 
CONDITIONS 


CALL  SUBROUTINE  EXTREM 


t 


PRINT  THE  FINAL 
PARAMETERS  AND 
RESIDUALS 


Figure  B2.   MAIN  Program  Functional  Flow  Chart 


57 


START 


READ  AIRCRAFT 
PARAMETERS 


I 


READ  AERODYNAMIC 
COEFFICIENTS 


PRINT  AERODYNAMIC 
COEFFICIENTS 


RETURN 


Figure  B3.   Subroutine  DATAIN  Function  Flow  Chart 
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START 


ASSIGN  VALUES  TO  COUNTERS  AND  FLAGS 


CALL  SUBROUTINE  FN 


INCREASE  STAGE  COUNTER 


COMPUTE  FIVE  ORTHOGONAL 
SECONDARY  DIRECTIONS  BY 
GRAM-SCHMIDT  PROCESS 


<^  ITERATION  DO  LOO? 


CALL  SUBROUTINE  FN  TO 

DETERMINE  CRITERIA  FN 

VALUE  AT  X-DX 


A 


T 


"INITIAL  ARGUMENTS 
OUTSIDE  BOUNDARIES" 


RETURN 


Figure  B4 .   Subroutine  EXTREM  Functional  Flow  Chart 
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CALL  SUBROUTINE  FN  TO 
DETERMINE  CRITERIA 
FUNCTION  AT  X  +  DX 


REDUCE  STEP 
SIZE 


FIT  PARABOLA 


CALL  SUBROUTINE  FN  TO 
DETERMINE  APPROXIMATE 
MINIMUM    (FOPT) 


REDUCE  STEP 
SIZE 


SET  NEXT  ORTHOGONAL  SEARCH 
DIRECTION  AND  STEP  SIZE 


DETERMINE  NEW  MAIN  DIRECTION 
AND  STEP  SIZE  FOR  THE  NEXT 
STAGE  WITH  X  =  XFOPT 


Figure  B4 .   (Continued) 
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START 


SET  SEARCH  VARIABLES 


CALL  SUBROUTINE  CONST 


LI=LI+1 


T 


CALL  SUBROUTINE  CLOKUP 


CALL  SUBROUTINE  SSEQNS 


CALCULATE  VALUE  OF  CRITERIA  FUNCTION 
&  INCREASE  ITERATION  COUNTER  BY  ONE 


RETURN 


Figure  B5.   Subroutine  FN  Functional  Flow  Chart 
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START 


CALCULATE  SPIN 
EQN.  VARIABLES 


RETURN 


Figure  B6 .   Subroutine  CLOKUP  Functional  Flow  Chart 


START 


DETERMINE  INTERPOLATION 
INDICES  FROM  ALPHA  &  BETA 


INTERPOLATE  ALPHA  &  BETA 
DEPENDENT  COEFFICIENTS 


INTERPOLATE  ALPHA 
DEPENDENT  COEFFICIENTS 


RETURN 


Figure  B7.   Subroutine  CLOKUP  Functional  Flow  Chart 
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START 


i 


CALCULATE  THE 
RESIDUALS  OF  THE 
SPIN  EQUATIONS 


RETURN 


Figure  B8.   Subroutine  SSEQNS  Functional  Flow  Chart. 

D.   PREPARATION  OF  THE  DATA  DECK 

The  following  is  a  guide  for  preparation  of  the  data 
deck  used  by  computer  program  SPIN. 

CARDS  1  and  2  -  The  initial  conditions  on  GAMDOT  (rad/ 
sec),  R(ft),  RDOT (ft/sec) ,  PHI (deg) ,  THETA(deg),  PSI (deg) , 
and  DELTAA(deg)  on  card  one  and  DELTAR(deg),  DELTAE(deg), 
THRUST (lbs)  and  ZDOT (ft/sec)  on  the  second.   Data  is  input 
as  floating  point  according  to  the  format  (7F11 . 5/4F11 . 5) . 

CARD  3  —  The  controlling  parameters  for  Subroutine 
EXTREM.   DFMAX,  DXMAX,  FMAX,  LMAX  and  IW  are  read  in  ac- 
cording to  the  format  (  3E10 . 2 , 15 , 12)  . 

CARD  4  -  Aircraft  identification  name  (e.g.,  CONFIGURA- 
TION A)  starting  in  column  one. 

CARDS  5  and  6  —  The  mass  and  physical  characteristics 
of  the  aircraft,  atmospheric  density  and  gravitational 
acceleration  values.   Card  5  contains:   SSAM  (mass,  slugs), 
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2  — 

S    (wing    area,    f t    ) ,    B    (wing    span,    ft),    CBAR    (c,    ft)    and    G 

2 

(gravitational  acceleration,  ft/sec  ).   Card  6  contains: 

2  2  2 

XI  (I  ,  slugs-ft  ) ,  YI  (I  ,  slugs-ft  ) ,  ZI  (I  ,  slugs-ft  ), 
x  y  z 

2 
ZXI  (I   ,  slugs— ft  )  and  RHO  (atmospheric  density,  slugs/ 

3 

ft  ) .   Data  on  each  card  are  input  according  to  format 

(5F16.8)  . 

CARDS  7—412  —  (alternating  three  card  series) .   These 
cards  contain  the  tabulated  aerodynamic  coefficients  that 
are  dependent  on  both  angle  of  attack  and  sideslip.   The 
first  card  in  the  series  is  as  follows:   in  columns  1&2 
the  coefficient  number  (I),  columns  3&4  blank,  column  5 
the  index  on  beta  (IB) ,  columns  6—10  blank  and  the  remain- 
ing 70  columns,  in  10  column  increments,  are  the  coeffi- 
cients indexed  on  alpha.   The  second  card  of  the  series  is 
divided  into  eight  10  column  increments  and  the  third  card 
is  divided  into  four  10  column  increments  of  aerodynamic 
coefficients  indexed  on  alpha.   Alpha  is  indexed  19  times 
for  each  beta  index  which  corresponds  to  a  range  of  alpha 
from  0  to  +90  degrees,  in  5  degree  intervals.   Beta  is 
indexed  nine  times  which  corresponds  to  a  range  of  beta 
from  —40  to  +40  degrees,  in  10  degree  intervals. 

The  first  set  of  three  cards,  within  a  particular  coef- 
ficient grouping,  corresponds  to  an  angle  of  sideslip  equal 
to  —40  degrees,  where  the  first  coefficient  in  columns 
11—20  corresponds  to  an  angle  of  attack  equal  to  0  degrees, 
columns  21—30  corresponds  to  an  angle  of  attack  equal  to 
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5  degrees  and  so  on  in  5  degree  increments  of  alpha  until 
90  degrees  is  reached  in  columns  41—50  of  card  three. 
The  second  set  of  three  cards  contains  coefficients  over 
the  alpha  range  for  beta  equal  to  —30  degrees.   The  remain- 
ing series  have  coefficients  for  beta  equal  to  — 20°,  —10°, 
0°,  +10°,  +20°,  +30°  and  +40°.   The  next  grouping  of  27 
cards  (nine  three  card  series)  contains  another  aero- 
dynamic coefficient  over  the  range  of  alpha  and  beta  and 
so  on  in  27  card  groupings  until  all  15  alpha  and  beta 
dependent  coefficients  have  been  prepared.   Format  for  the 
three  card  series  is  (12 , 2X, II , 5X, 7F10 . 7/8F10 . 7/4F10 . 7) . 

The  input  sequencing  of  groupings  is  very  important  as 
the  machine  computation  is  dependent  upon  the  assignment 
of  a  particular  tier  in  the  three— dimensional  COEFAB  array 
to  the  proper  aerodynamic  coefficient.   Table  Bl  gives  the 
proper  assignment. 

CARDS  413—429  —  (alternating  three  card  series) .   These 
cards  contain  the  tabulated  aerodynamic  coefficients  that 
are  dependent  on  alpha  only.   The  first  card  in  the  series 
is  as  follows:   columns  1&2  the  coefficient  number  (I), 
columns  3—10  blank  and  the  remaining  70  columns,  in  10  col- 
umn increments,  are  the  coefficients  indexed  on  alpha.  The 
second  and  third  cards  of  the  series  are  the  same  as  de- 
scribed above  for  cards  6—411.   Again  alpha  is  indexed  19 
times  corresponding  to  a  range  of  alpha  from  0  to  90°  in 
5°  increments.   Columns  11—20  of  card  one  corresponds  to 
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the  coefficient  for  alpha  equal  to  zero  degrees,  columns 
21—30  for  alpha  equal  to  5°,  and  so  on  in  5°  increments 
until  alpha  equal  to  90°  is  reached  in  columns  41—50  of 
card  three.   Format  for  these  three  card  series  is  (12, 8X, 
7F10. 7/8F10. 7/4F10. 7) .   As  before,  the  input  sequence  of 
the  coefficient  groupings  is  very  important  for  proper 
interpretation  by  the  program.   The  assignments  are  given 
in  Table  Bl. 
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TABLE  Bl.   Coefficient  Assignment 


Storage 

Coefficient 

Interpolated 

Coefficient 

Location 

Print  Out 

Location 

Cl 

COEFAB  (  1,  ,  ) 

1 

CR(  1) 

C 
m 

COEFAB  (  2,_,_) 

2 

CR(  2) 

C 

n 

COEFAB  (  3,  ,  ) 

3 

CR(  3) 

C 
Y 

COEFAB  (  4,  ,  ) 

4 

CR(  4) 

C 

X 

COEFAB  (  5,  ,  ) 

5 

CR(  5) 

C 
z 

COEFAB  (  6,  ,  ) 

6 

CR(  6) 

Cy6r 

COEFAB  (  7,_,_) 

7 

CR(  7) 

C,  x 
l<5r 

COEFAB  (  8,  ,  ) 

8 

CR(  8) 

n6r 

COEFAB  (  9,  ,  ) 

9 

CR(  9) 

y6a 

COEFAB  (10,  ,  ) 

10 

CR(10) 

r 
16a 

COEFAB  (11,  ,  ) 

11 

CR(ll) 

n6a 

COEFAB  (12,  ,  ) 

12 

CR(12) 

z6e 

COEFAB  (13,  ,  ) 

13 

CR(13) 

C  x 

m<5e 

COEFAB  (14,  ,  ) 

14 

CR(14) 

C   X 

x6e 

COEFAB  (15,  ,  ) 

15 

CR(15) 

C 

YP 

COEFA   (  1,  ) 

16 

CR(16) 

Cl 

COEFA   (  2,  ) 

17 

CR(17) 

P 

C 
n 

P 

COEFA   (  3,  ) 

18 

CR(18) 

C 

COEFA   (  4 ,  ) 

19 

CR(19) 

Cl 
r 

COEFA   (  5,  ) 

20 

CR(20) 

C 

n 
r 

COEFA   (  6,  ) 

21 

CR(21) 

C 

COEFA   (  7,  ) 

22 

CR(22) 

q 

c 

m 

q 

COEFA   (  8,  ) 

23 

CR(23) 

c 

X 

q 

COEFA   (  9,  ) 

24 

CR(24) 
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E.       PROGRAM    LISTING 

£  $  #  $  V  $  :fr  :<.  $  ^-  £*:£££#  y  :;::*:;■:*###£:<;#  # 

c  * 

C  MAIN  PROGRAM  * 

C  * 

C  PURPOSE     SERVES  TO  SET  THE  INDEPENDENT  VARIABLES  TC  * 

C               BE  SEARChEC  OVER  ANC  CONTROLS  ThE  REACING  CF  * 

C               THE  INITIAL  CONDITIONS  ANC  OPTIMIZATION  * 

C               CONTROL  PARAMETERS.   ALL  STRUCTURE  FOR  THE  * 

C               PRINTED  OUTPUT,  rtilTH  THE  EXCEPTION  CF  THE  * 

C              OPTIMIZATION  RESULTS,  COMES  FROM  THIS  * 

C               PRCGRAM.  * 

C  * 
Q $*  :M  4* 4*$ $$$$$$ ^ $£$$$$£$ $4*4$* £$$*$$ $$*$# $$*$##$ *$#$$$$$$$£ 

IMPLICIT  REAL*S( A-F,0-Z) 

EXTERNAL  FN 

DIM  ENS  ION  Xlll),DX(ii),A(lif14) 

COMMON/ SP I N/G AMD GT ,GAMMA,RDCT,Phl,THETA,PSI,CELTAAfDEL 
1TAE,DELTAR,THRC;ST,V,COEFaB(15,20,10)»COEFA(20,9),SSAM, 
2S,6,CdAR,XI,YI,  ZI  ,ZXI,&,RHO  ,All,AU,A13,A2l,A22,A23,A3 
=  i,A32,A33,VSw,VXl,VYl,VZl,P,C,RZ,PEAR,QbAF-,F-dAR,HRhCS, 
ACR(24),RES(6i  ,  ALPHAD , ALPHA R ftiET AC , &ETAR , SEE TA ,C8ETA ,SA 
5LPhA,CALPHA,SGAMMA,CGAMMA, SPSI, C P S I , SThE T A , CTriETA , SPH I 
6,CPHI ,R,XFCRCE,XFORl , Y FORCE , YFCR 1 , ZFORCE  ,  ZFCR1 , XMOM , YM 
70M,  ZM0M,ZDQT,NPRINT,NAME(4) 
C 

C  DATA  INITIAL  STEP  SIZES  OF  THE  INDEPENCENT 

C  VARIABLES 

C 

DATA  OGAM/0.1/,DPSI/.01/,DZDQT/5 .0/ , CPH I / .C 1/ , OR/ 5 . C/, 
lDTHETA/.31/,DDELE/C.l/,DTHRLS/lCC.C/,DRDCT/5.0/ 
C 

C  READ  THE  INITIAL  CCNDITICNS  CN  THE  INCEFENCENT 

C  VARIABLES  AND  OPTIMIZATION  CCNTROL  PARANETERS 

C 

RE AD (5, LOO)  GAMDCT , R, ROOT ,PH I , TF ETA , PS  I ,CELTAATDELTAR, 
1DELTAE, THRUST, ZDDT 

REAC13, 101 i  DPHAX,CXMAX,FMAX,LMAX,  lb 
C 

C  SET  PRINT  INDICATOR  FOR  THE  AERODYNAMIC 

C  COEFFICIENTS:       NPRINT=0,  SUPPRESS  PRINTING 

C  NPR1NT  =  1,  PRIM  COEFFICIENTS 

C 

NPRINT=0 
C 

C  SET  THE  NUMBER  OF  INDEPENDENT  VARIABLFS  (K=)  ANC 

C  ASSIGN  EACH  NAME  TC  THE  X  ARRAY  WITH  THE  CORRES- 

C  PCNDING  STEP  SIZE  TO  THE  DX  ARRAY. 

C 

K  =  6 

X( 1)=GAMD0T 

X(2)=R 

X(3)=ZDCT 

X(4)=PSI 

X(5)=THETA 

X(6)=PHI 

DX<  1)=DGAM 

CX<2)=DR 

CM3)=DZDCT 

CX(4)=0PSI 

DX(5)=DTHETA 

DX6)=DPHI 
C 

C  CALL  TC  SUBROUTINE  CATAIN  TC  READ  MASS,  PHYSICAL 

C  CHARACTERISTICS  AND  AERCCYNAM1C  COEFFICIENTS 

C  OF  ThE  AIRCRAFT. 

C 

CALL  DATAIN 
C 

C  WRITE    OUT    THE     INITIAL    CONDITIONS 

C 

WRITE (6, LLC)     NAME 

kFITE(6,102) 
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*  MAIN  PROGRAM  CONTINUED  * 

kRITE(6,103) 

WPI7E(5,1U4)  GAMDCT,R,RDOT,PHl,TFETA 

WPIT6(o,105) 

kRU£(  6,104)     PSI,DELTAA,DELTAR,DELTAE,THPLST 

WRIT£(6,  106) 

kRITE(6,119)  200T,  PHC  ,  DEM  AX  ,  CXM  AX  ,  FMAX 

kRITE( 6,  107) 

kRITE(6,  103)  LMAX,  Ik 

kPITEU  ,109) 

CALL  THE  OPTIMIZATION  ALGORITHM  TO  FINE  THE 
STEADY  SPIN  PARAMETERS. 

CALL  EXTREM(FN,K,X,0X,A,DFMAX,DXMAX,FMAX,LMAX,FOPT,Ik) 

WRITE  CUT  THE  FINAL  STEADY  SPIN  PARAMETERS. 


k 
k 

k 

k 
k 

v. 
k 
h 
k 

D 
k 
k 
S 
F 
F 
F 
F 
1» 
2« 

104  F 

105  F 
1« 
2« 

1C6  F 
1' 
F 


ICC 
101 
102 
102 


107 
ICE 
109 
11C 
1  11 
112 


F 

F 

P 

F 

F 

1« 

2' 

112  F 

IT 

2« 

114  F 
1* 
2' 

115  F 
1« 

F 

F 

IT 

IF 

F 

1' 

119  F 

12C  F 

E 


116 
117 


118 


PITE(6 

RITE(6 

P1TE(6 

PITE(6 

RITE( 6 

RITE(6 

RITE(6 

PITE(6 

RITE(o 

RITE( 6 

RITE(6 

C  2  1  = 

PITE(6 

PITE(6 

TCP 

CFMAT 

CFMAT 

CRMAT 

CRMAT 

FPi1  ,T 

(FT/SE 

CRN  AT 

CRMAT 

DELI  AT: 

(DEG) ' 

CRMAT 

C>MAX« 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

TFETA' 

(RAD/S 

CRMAT 

1C3,  'k 

(LPS)  • 

CRMAT 

F'  ,/T5 

(FT/SE 

CFMAT 

fT67,  • 

CFMAT ( 

CRMAT 

6  9  ,  '  Z  F 

15.5/  , 

CRMAT 

CF  THE 

CRMAT 

CRMAT 

NC 


,110)  NAME 

,111) 

,112) 

,104)  ALPHAC,6ETAD,GAMD0T,TFETA,PFI 

,112) 

,104)  PSI, GAMMA, R, THRUST, RFC 

,114) 

,104)     V,  VX1,VY  l,V/ZlfP 

,115) 

,116)     Q,RZ,RDCT 

,117)     XFORCE, YFCRCE,ZFORCE,XMQM,YMCN  ,ZMCM 

1,6 

,120)     I,RES(I) 

,11£)    FOPT 


(7F1 
(3E1 
(T67 
(//, 
102, 
C)«  , 
(T54 
(//, 
»,T1 
tT91 
(//t 
i  TLC 
(//, 
(T57 
(51/ 
(•  1' 
(T67 
(//, 
,T1C 
EC)  • 
(//, 
HC  •  , 
,/) 
(//, 
E,«  ( 
C)«  , 
(//, 
(RAD 
T54, 
(/// 
CRCE 
T71, 
(/// 
RES 
(T3B 
(//, 


1.5/4F1 
0.2, 15, 
,'-  -  - 

T5&, »GA 

■THETA* 

i 


1.5) 
12) 

INI 
MDGT 
f/T5 


EG) 
) 

I*  ,T 
1ST1 
•tT  L 
156,' ZDCT' , 
«,/T 
AX'  , 
,13, 


TIAL  CCND 
T71,fR« 

(RAD/S 


ITICNS      -') 

191, 
CG)  '  ,T7 


,5F11.5 

T58, ' PS 
02, ' THR 
, • (DEG) 


,T£C, 'ROOT  '  , 
EC)  '  ,To9,  * 
«  ,T102  ,« (DEC)'  ,/) 


(D 


T55,  • LM 
,  I 5,T6£ 


)  ,T67,  • 
,8(/)  ,T 
,  • 

T58, • AL 
2,  'PHI  • 
tTSL, '( 
T59, «  PS 
/T56, « ( 

T61,*  V 
FT/SEC) 
T1C1, « ( 
T61,»  Q« 
/SEC) 
^Fl  i  . 
,T69, 
=  «,FU 
• ZMOM= « 
,T58, 'S 
IDUALS1 
,2F11  .5 
T69,  'RE 


72,4 
FIN 
PFAD 
,/T5 
DEG) 
I"  fT 
DEG) 


69,  • 
,/T5 
02,  ■ 
T69, 
55,' 
T69, 
/) 

-  EX 
A4,/ 
AL  P 
•|T6 
8,'  ( 
1  ,Tl 
69,  ■ 
'  ,T6 


DELTA 
7f  •  (D 
(LBS) 

•  R  H  Q  < 
(FT/S 

•  Ik1  , 


A'  ,T80, ,DELTAR,fT91, 
EG)  •  ,T69, «  (DEG)'  ,T80 
«  ,/) 

jTSC'DFMAXSTgi, 
EC)  '  ,/) 
/) 


TREM     ITERATIONS    - 

/  ) 

ARAME 

9,«  BE 
DEG)  • 
0  1 ,  •  (  C  E  G  )  '  ,  /  ) 

GAMMA 


_    _  i 


,//) 


TERS 

TAC',T80, 


-•  ) 


T69, ■  (DEG)  « 


MDCT' ,T91, 
,T78, 


C,  ■THRUST1, 
(DEG)  ' ,T&0,  '(FT)'  ,T91, 


,T70  ,  ' VX1  ■  ,T£ 
•  ,T68,'  (FT/SE 
RAD/SEC) «  ,/) 


1,  «VY1«, 
C)  '  ,T79, 


T91,  •VZl'  ,T1C5, 
•  (FT/SEC)  •  ,T69, 


T70,,RZ«  ,T8  1  ,  'ROOT' ,/T57,«  (RAD/SEC) 


5, 
«X 


,/) 


T79, ' (FT/SEC) 

5(/  )  ) 

FCRCE=' ,F15.5/,T69,'YFCPC 

5/,T7  1, »XMOM= 

,F15.5,2(/  )) 

UM    CF    THE    SWU 

,E17.5) 

,2X,3E11 .4) 


•  ,F15.5/,T71 
ARES' ,/T58, 


E=«,F15.5/, 

,• YMCM=«  , 


S(  '  ,  I  1,"  ) =■  ,E 


I 


) 


69 


£**:***************************#********************$******** 

c  * 

C  SL6PCLTINE  CATAIN                                           * 

C  * 

C  PLRPCSE     READS  THE  MASS  AND  CIMENSICNAL  CHARACTERISTICS 

C  OF  THE  AIRCRAFT  AS  WELL  AS  THE  AEPCCYNAMIC    * 

C  COEFFICIENTS.                                  * 

C  * 

^^^4^A^^^*^^^*V;«^:'/.^^i*^**^^^^^*i^**A;^*^^A^**^^H--************^ 


c 

C 
C 

c 

c 


c 

c 

C 

c 


SL6RCUTINE  DATAIN 

I  PPL IC IT  RBAL-8( A-P,3-Z) 

CCf'NCN/SPIN/GAMDCT,GAf-TA,RDOT,PP  I  ,  TPETA,  PS  I  ,CELTAA,CEL 
IT  A E, DEL  TAR, IHRUST, V,CQEFAB(  15 
2S,&,C6AR,XI,YI,ZI, ZXI,G,RHO,A 
3 1,  A32 » A33 » VSQ, VX 1, VY 1, V Zl ,P,G 

4CR  (24)  1  RES (6)  ,  ALP  HAD,  ALPHAR,cETAL  ,  GETAR  ,  S  E  E  TA  ,CbE  T  A  ,  SA 
5LFPA,CALPHA,SGAKMA,CGA^A,SPSI,CPSI,STHElA.CTHETA,SPHi 
6,CFHl,R,XFCi-CE,XFCRl,YrCRCE,YFGRl,ZF0RCE,ZFCRl,XFCN,Yi^ 
7CM,ZMOf',ZDCTf  NPRINT,NAME(4J 


,tCi  1C) ,COEFA(2C,9) , SSAN, 
11,A12,A 13,A21,A22,A23,AS 
,RZ,FEAR,wEAR,P6AR,FRFCS, 


READ  IN  THE  CCNF I  CURAT  I  ON  NAfE, 
PARAMETERS  OF  TFE  AIRPLANE,  AIR 
GRAVITATIONAL  ACCLEPATICN  VALUE 


MASS  ANC  PHYSICAL 
DENSITY  AND 


READ<5,100)  NAME 
REAC(5,1C1)  SSAM, 
READ(5,  101)  XI, YI 


S,E,CBAR,G 
,ZI  ,ZXl  ,RHO 


REAL  TFE  AERODYNAMIC  COEFFICIENTS  THAT  ARE  A 
FUNCTION  OF  ANGLE-CF-ATT ACK  AND  SIDESLIP 

DC  1  J =1,1 35 

READ(5,1C2)  I  ,  IB,  (CCEFABU  ,  IA,  IB)  ,  IA  =  1,  19) 


READ  TFE  AERODYNAMIC  COEFFICIENTS 
FUNCTION  OF    ANGLE-QF-ATTACK  CNLY 


THAT  ARE  A 


It 
If 


DC  2  J  =  l,6 
READ( 5,  103) 
CC  3  J=l,3 
REAC(5f 1C4) 
CC  4  1=1,  1? 
CC  4  J=1,1C 
CCEFABt I,2C, J)=O.C 
CCEFABl I ,J  ,  LO)=O.C 
CC  5  K=l,9 
CCEFA( 2C,K)=0.0 


(COEFA(  IA 
(COEFA(  I  A 


I  ) 

I  ) 


IA  =  1 

IA  =  1 


IS) 
19) 


CHECK  PRINT  INC1CAT0R  FOR  ECFG  CHECK  CF  DATA 


IF< 

CC 

hRI 

WR] 

hPI 

kRI 

WRI 

hRI 

N  =  C 

CC 

hRI 

N  =  N 

CCN 

N  =  C 

WRI 

VhP  I 

WRI 

DC 

WFI 

N  =  N 


NPR INT. EG 
7  1=1,15 
TE(6 


TE(6 
TE(fc 
TElc 
TE(6 
TE(6 


105) 
11C) 
106) 

107) 
108) 
1C9) 


0)  GO  TO  9 
NAME 


6  IA=1,19 

TE(6,LU) 

+  5 

T  INUE 

TE.(6,  10  5) 
TE(6  ,113) 
Tt(c,U2) 
6  1=1,  19 
TE(6,111) 
+  5 


N, (CCEFABt  I, IA,  IB) 


NAME 


N,(CCEFA(I,J),J=1,9) 
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-'-  O,  -X-  »'- 


#^: 


:*^Jt^*4);:/c)lt: 


££**$£ 


1C 
10 
10 
1C 
1C 
1C 
10 
10 


10 
10 

11 
11 

li 


11 


9  P 


F 

F 

F 

F 

F 

F 

F 

F 

1" 

2' 

8  F 


F 

F 

F 

F 

1« 

21 

3' 

2  F 

E 


E-TURN 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

CR^AT 

CRMAT 

CRMAT 

CRMAT 

-10'  , 

+  4C  ) 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

CRMAT 

CCEF 

92, «C 

(CEGR 

CRMAT 

t\C 


SUBROUTINE  CATAIN  CONTINUED  * 


(4A4) 
(5F16 

(12 
(12 
(12 
('  1 
(T6 
(T29, 


>  * 
,  1 

>3 , 


.8) 

X, I 1,5X,7F1C. 

X,7F10.7/8F10 

X,7F11.7/7F11 
9(/)  ,T£2,4A4) 

•  COEFFICIENT 
•BETA  •  ,T40i  ■- 


7/8F10.7/4F10. 
•7/4F10.7) 

.7/5F  11.7) 


•  ,12 

40  « 


//) 


T77,  '0'  ,T35,  «  +  10* fT94 


,  T  4  9  , 

•+2C 


i  _ 


30*  , 
,T103, 


7) 


T5£, 

'  420 


(T2 
(T2 
(/) 

(T3 
(T2 

L8'  , 

OEF 

EES) 
(/) 


7, ! (DEGREES) ■ ) 
9,«  £LPPA'  ) 


C,  I2,4X,9F9.5,/ ) 

9,'  ALPHA*  ,T38  ,'CCEF  1 6 •  ,T47 , • CC E F 
19*  ,T7A  ,  'CCEF  20'  ,T32 


T65, 'CuEF 
22'  , T101  , 
',/) 


•-20' ,T67 
«,T112, 


17' ,T56, 
21'  , 


CCEF  23  •  , ri  1C,  'CutF  24',/T27, 


71 


J**.**********************,********************************* 

* 
SLBRCUTINE  EXTREM  * 


PLRPCSE 


C** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C** ********* 


THE  G 

MAIN 

FOB  A 

SUBRu 

OP  SE 

ARCUM 

THE  M 

A  PAR 

VALUE 

FN. 

DETER 

PROCE 

SIZES 


PTI  MI 

PR-JGK 
PART 
UTINE 
ARCH 
ENT  L 
IN  I  Mb 
ABOLA 
POIN 
SEC  ON 
MINED 
SS  AN 
t  UNT 


ZATI 
AM  T 
ICUL 

IS 
BY  T 
1ST 
H  AL 

ThR 
TS  D 
DARY 

BY 
C  CO 
IL  A 


CN  A 

HAT 
AR  S 
PROV 
HE  U 
AND 
CNG 
GUGH 
ETER 

DIR 
A  GR 
NT  IN 

M  I N 


LGOR  I 

RETUR 
TEADY 
IDED 
SER  T 
THE  A 
THIS 

THRE 
MIKED 
ECTIC 
AM-SC 
UES, 
I  MUM 


TPM  C 

NS  TH 
SPIN 
WITH 
HRCUG 
LGCRI 
LINE 
E  C^  I 
BY  C 
NS  GF 
FMICT 
V*  I  Th 
IS  FG 


ALLEC  FPCN  THE    * 

E  MINI  MM  FCUND   * 

SEARCH.   THE     * 

AN  INITIAL  LINE   * 

h  ThE  CALLING     * 

THM  C El  ERMINES    * 

BY  INTERFCLATING  * 

TERIA  F INC T ION    * 

ALLS  TC  SLdROUTINE 

SEARCH  ARE        * 

ORTPCGCNALIZATION 

VARYING  STEP      * 

UNO  .  * 


•  -r  *r  -r  t*  *y*  i* 


:^^i 


.  u.  ~±   O 


**********:?***** 


■v  -r  *f  • 


;  *.  ir  i.:  if  if.  if.  *  *  &  4  £  if  jj: 


c 
c 

C 

c 
c 


c 

c 

c 


c 
c 

c 
c 


SLERCLIT  INE  EXTREM  (  F  ,K  ,  X  ,DX  ,  S  ,  DFM  AX  ,  CXMAX  ,  FM  AX  ,  LMAX  » 
1FCPT, IK) 
IMPLICIT  REAl_*S(  A-P,Q-Z) 
CIMENSICN  X(l) ,DX(1) , S ( 1 1 ,1) 

ASSICN  COUNTER  AND  FLAG  VALUES 

L  =  C 

L  1  =  1 

N  =  C 

DC  1  1=1, K 

S  ( 1 , 1 )  =  X  ( I  I 

1  S(I,2)=X(IJ-DX(I) 

CALL  "IC  SU3R3UTINE  FN  TO  CHECK  IF  THE  INITIAL 
CONDITIONS  EXCEED  THE  BOUNDARIES  OF  ANGLE-CF- 
ATTACK  AND  SIDESLIP 

CALL  F( X,F2,LI ,N) 

FE  =  F2 

I  F  (  L  I  .GT.l )WRITE(6,2) 

2  FCPMATI  IX,  • INITIAL  ARGUMENTS  OUTSICE  BCUNCAFIES') 
2  I  F(KC  .GE.K.CR.KC .LT.O.OR.L . EQ.G)KC=0 

KC=KC+I 
S(  1 ,3 J  =  C. 
CC  4   1  =  1, K 
S(I  ,4)=S(I  ,1)-S( I ,2) 
4  S(i,3)=S(l,3)+S(I ,4)**2 
S(  1,3)=CSCRT(S<1, 3) ) 

IFUABS(IW).GE.3)WRITE(6,23 J L,N ,S ( 1, 3 ) , F2 , 
1(  I,X(I) »I,DX( I) ,I  =  i,K) 

CHECK  TEE  STOPPING  CONDITIONS 

IF(L.GE.IA6S(LMaX) .OR.S(l, 3) .LT . CXMAX .OR . CA ES ( FF-FOPT) 

l.LT.DFMAX.AND.L.GT.O.CR.(LI.GT.l.Ai\D.L.Ev».C).CR.FOPT. 
2LE  .FMAA  .ANC.L.GT .0 JCO  TO  22 
IF(K.E^.i}GCTC9 

COMPUTE  THE  ORTHOGONAL  SECONCARY  DIRECTICNS  BY  A 
GRAM-SChMICT  CRTHOGON AL IZAT I CN  PROCESS 

CC  8   J  =  2,K 

KC=-2  +  J-*KC 

IP(KD.GT.K)KD=KD-K 

S(J,3)=0. 

DC  7   1=1, K 

S(I,J+3)=0. 

IF ( I .EC  .KC)^{ I, J  +  3)  =  S(  1,3) 

JM=J-1 

CC  6   JK=1,JM 
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c 
c 
c 


e 
9 

1C 


SUBROUTINE  EXTREM  CONTINUED 

6  S( I,J  +  3  )  =  S(  I ,J+3)-S(KC, JK+3)*S(  l,3)/S(JK,3)*S(I,JK+3) 
1/S( JK,3J 

7  S(J,3)=S( J,3)+S(  I  ,J+3)**2 
S(  J,3)=CSQPT(S(J,3 ) ) 
IF(SC  J,  3)  .LT.1.D-3OG0TC3 
CCNTINUE 
DC  10   1=1 f K 
S( I , 2  J  =  S( I  ,1) 

L  =  L  +  1 
FF=FCPT 

STAGE  CO  LOOP 

DC  2  1  M=1,K 
DC  1  1  1  =  1, K 

11  S  (I  ,M+3)=S(  1 1  .v  +  3)/S  (M,3)*DX (M) 
IF(IW.GT.0)LI=3 

12  IF  (IW.GT.O  )L I  =  L 1-1 
LJ  =  LI 

DC  13  1=1, K 
X(  I)  =  S(  I,  U-S(I,M  +  3) 
S(I  ,M+3)=S(I,1)-X( I  J 

CALL  SUBROUTINE  FN  TO  DETERMINE  THE  CRITERIA 
FUNCTION  VALUE  AT  X-DX 


1 


CALL  F(X,F1,LI,N) 
BC  =  1  • 

14  C  C  15  1  =  1, K 

X(  I)=S(  1,1 )+S(I ,N  +  3)/BC 

15  SCI ,M+3)=X( i)-S( I  ,  1) 

IF  (DABS  (BO  J  .GT.l.  i)GQT.020 

CALL  SUBROUTINE  FN  TO  DETERMINE  THE  CRITERIA 
FUNCTION  VALUE  AT  X+DX 

CALL  F(X,F3,LJ,N) 

IF  CUT^IDE  THE  BOUNDARIES  GC  EACK  AND  ACJLST 
THE  STEP  SIZE 

IFCLI+LJ.EC.4  )G0TC12 
IF  (LJ .CT.2 )BC=-4  . 
IF(LI  .GT.2)B0=+4. 
IF(LI.GT.2.0R.LJ.CT.2 )GCT014 

16  ST=C. 

INTERPOLATE  A  PARABCLA  TFROUGF  THE  THREE 
FUNCTIONAL  VALUE  POINTS:  X-DX,  X,  X+CX 

CC  18  1=1, K 
X( 1)=S(  1,1  ) 

IF(DABS(S(I,M+3)).LT. UD-30)G0TC1£ 
S ( I ,M+3 )=S(I , K+3)/LI 

1F(CABS(2.*F2-F1-F3) . LT. 1.0-30) CCTC18 
X(  I)  =  S(  1,  L)+S(I,M+3)/DABS(F  1-2.  *F2+F3 )  *  (  F3-F 1 )  / 
1IS1GN(2  ,LNAX) 

18  ST  =  ST  +  (  M  I)-S(  1,1)  )**2 

IF(16.*ST.LT.DX(M)+*2 )DX ( M ) =DX ( M  )  /4 

1F(ST.LT.4CC.*DX(M)**2  .AND.  DABS  (2.*F2-Fl-F3).GE.l.D-30 

.1)GC  TO  2C 
rr  ] g  y  —  i  v 

IF (DABS (S( I,M+3) ) .LT. 1 . C-30 ) GOTC 19 

INTERPCLATE  FOR  THE  APPROXIMATE  MINIMIS  POINT 

X(  I)=S(  I,i)+DSIGN(S(I ,M+3)  ,  ( F 2-F  1 ) /S( I , M  +  2 ) ) * 
1  ISIC.\(2C,LMAX) 

19  CCNTINUE 
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4* 


c 

c 

c 


*  $  $  3^  if.  ;J  y:  £  £  >lc  *  3i  i-  3^  *  :!:  *  #  ^  £  *  *  aft.*  *  *  sjt  ajt  *  *  *  *  *  *  *  *  *  *  *  *  **  *  *  *  *  *  *  **  *  *  *  *  *  #  *  i- 

SUERCUTINE     EXTREM    CCNTINUEC  * 

CX  (N)=DX(M)*2. 
2C    LI-+1 

EO-6G 

IF (CABS (EG)  .GT .1 .1 )DX(M)  =  DX(K)/2  . 

DETERMINE    THE    MNIMUM    VALUE 

CALL    FO,FCPT,LI,N) 
IF(LI.G1  .1  )LI=10 

SET    THE    NEXT    CPTHCGONAL    SEARCP    DIRECT  ICh    £ND 
STEP    SIZE 

IFUSIGMl  ,LMAX)*(F0PT-F2)  .LT.-C£ES(FE-F2)*4.  .AND. LI  . 
1NE.10)LI=2 
IFtLI.GT.l .AND.DABS(BO) .GT . 1 . 1 ) GGT014 
IFtLI .GT.l ) GOTO! 6 

FE  =  F2 

F2=FOPT 

DC    2  1    1=1, K 

21  S(  I, L)  =  X(I  ) 
GCTC3 

22  IFdABS  (IiO  .EC.2)  WRITE(6,23)L,N,S(  1,3)  ,  F2  , 
1 ( I,X( I  )  ,1 , CX(  II , 1=  1  tK  ) 

23  FCRMAT(//,T57, 'STAGE    NO  . ■  , I  5 ,20X  ,  ' TR  I AL    NC . • , I  6, //T 57, 
1,CL--',E15.5,//T57,  'FUNCTION    VALLE=',E15.5,///(T57, 
2,£R(,,I2,,)=I,E15.5,5X,,DS(',I2,,)  =  ,,E15.5H 

RETURN 
ENL 
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c** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c** 


A  %  $  ^  :,'  v  «  $  *  *  *  £  *  *  :'-.  $  A  *  *  *  A  *  *  *  *  V  *  *  *  *  *  *  *  *  *  #  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


SLbRCUTINE  FN 


F I R  F  C  S  E 


CAL 
CF 

FN 

AND 

VAL 

ANC 

THA 

LIN 

IS 

COR 

LYI 


LED  BY  SLE'FOUTINE  EXTREM  FOR  C ET E RN  I  NAT  ICN * 
THE  CRITERIA  FACTION  VALUE.  SIERCUTINE  * 
MAKES  CALLS  TO  SUBROUTINE  CONST,  CLOKUP,  * 
SSEQNS  TO  CETERMNE  TFE  CRITERIA  FUNCTICN* 
UE.  A  CHECK  TO  INSURE  THAT  THE  VALUES  CF  * 
LE-OF-ATT/5CK  AND  SILESLIP  IS  PROVIDED  SC  * 
T  A  SEARCH  bEYGND  ThE  CCEPFICIEfvT  TABLE  * 
ITS  IS  PROHIEITEC.  IN  SUCH  A  CASE,  CCNTPCL 
RETURNEE  TO  SUBROUTINE  EXTREM  SC  THAT  THE  * 
RECT  STEP  SIZE  CAN  BE  ACJUSTEC  TC  VALUES  - 
NG  WITHIN  THE  TABLE.  * 


*^*^**^«^^^*^***>!«^****^*****5?c4:^**4:****iie*^***5t*3{<^si.*5jt*3}ej{!*^* 


C 

c 
c 

c 
c 
c 


SL 
IF 
CI 

cc 

1TA 
2S, 

31, 
4CR 

5LF 
6,C 
7CF 
CA 
R  = 
lL 
PS 
TF 
FF 


ERO 
PLI 
P  EN 
MMG 
E,D 

e,c 

A32 
124 
FA, 
FFi 
i  ZM 
NEC 
X(2 
CT  = 
I=X 
ETA 
I=X 


UTIN 

CIT 

SICN 

M/SP 

ELTA 

BAR, 

,A3^ 

)  ,  RE 

CALP 

,R,X 

OM,Z 

T  =  X( 

) 

X(3) 

(4) 

=X(5 

(6  J 


E    FN(X, F, LI  ,N) 

REAL*6( A-H,C-Z) 

X(  I) 
IN/GAMLCT  ,GAMMA,RCOT,PHI,THETA,PS  I,CELTAA,CEL 
R,THRUST,V,COEFA3(15,20,1C),COEFM2C,9),SS£N, 

ii,AU,A13,A21,A22,A23,A3 
,RZ,PEAR,QBAR,RBAR,HRFCS, 
S(6},ALPHAD,ALFHAR,BETAB,BEfAR,SEETA,CBETA,SA 
PA,SGAMMA,CGAMMA,SPSI,CPSI,STHETA,CTHETA,SPHI 
FCRCE,XFCRi,YFCRCE,YFOPl,ZFCRCE,2FCRl,XMC!SYVi 
DCT,NPRINT,NAME(4) 

l  ; 


XI, YI,  ZI , ZXI ,G,RHO,A 
,VSQ,VX1 , VY1, VZ1  ,P,C 


CHECK  TO  SEE  IF  THE  60UNCARIES  ARE  VICLATEC. 
SO,  CONTROL  IS  RETURNED  TO  SUBROUTINE  EXTREM 
THE  STEP  SIZE  IS  ACJUSTEC  CF  NEW  ARGUMENTS 
DETERMINED  THAT  DC  NOT  VIOLATE  THE  BCUMARIES 


IF 

Inhere 


CALL  CCNST 

IF(ALPHAD.LT.0.0.0R.ALPHAD.GT.90.  .  CR .BE T AC  .  LT  .-40  .  .CR. 
1EETAD.GT.4G. )LI=LI+1 
I  F  (LI .GT.l JRETURN 

CUOKUP  MAKES  THE  INTERPCLATEC  COEFFICIENT  VALUES 
AVAILAELE  IN  CCMMCN  MEMCFY 

CALL  CLCKUP 

SSECNS  MAKES  THE  RESIDUALS  AVAILABLE  Ih  CCFMCN 

MEMORY 

CALL  SS EQNS 

CRITERIA  FUNCTICN  CETERMNATICN 

F=RES(  l)*RES( i)+RES(2)*RES(2)+RES(3)*RES(3)4RES(4)* 
1RES(4)+RES(5)*RES(51+RES(6)*RES(6) 

N  =  N  +  1 

RETURN 

ENC 
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*  *  *  $  *  ^  *  *c  $  v-  *  *<  *  *  *  *  *  *  *  -•*  *  *  *  *  £  ii  *  *  *  *  *  *  *  *  ^  *  *  *  *  *  *  4  4  *  *  £  *  * 

CLTINB  CCNST 


PURPOSE 


c 

C   SLBP 
C 

c 
c 
c 
c 
c 
c 
c  ■        * 

\^  51  -*  -y  3?  V  3$C  n  SfC  3$!  -,  *,<  -  -  -v  5J"  ^  -v  '.*  o-  ^  ^  V  v  5r  ^  ^  *r  ^v  -r  V  r  t  t  *^  »r  ->•  -"i"  ^r  V  1-  V  n*  *r  •**  *  3f  3p  »r  -V-  *?  V  t  n*  *  *  *  ^  '*■  T  '•* 


DETERMINES  VARICUS  EQUATION  PARAMETERS  NE 
IN  EVALUATING  THE  EOUATICNS  OF  MCTICN.  V 
OF  ANGLt-CF-ATTACK  AN'C  SILESLIP  ARE  PASSE 
TO  SUBROUTINE  FN  FOR  A  LIMT  CHECK  ML  AL 
TO  SUBROUTINE  CLOKUP  FOR  CCEFFICIENT  INTE 
PCLATION. 


***** 

* 

* 

* 

E  C  E  C  * 

ALLES 

C    * 

SC   * 

R-   * 


SIPRCU 
1  f-  F  L  1 C 
CCPNGN 
1TAE,DE 
2St 6,CB 
31.A22, 
4CR(24) 
5LPFA,C 
6,CPHI , 
7CNiZN!C 
7  HE TAR 
PFIR=P 
PSIR=P 
S1FETA 
CTFETA 
SPPI=D 
CFHI=D 
SFSI=D 
CFSI=D 
A12=-S 
A22=SP 
A33=CP 
A12=CT 
A22=SP 
A22=CP 
A11=CT 
A21=SP 
A21=CP 
HRFOS= 
VSC=<G 
V=DSQR 
F=A13* 
C=A23* 
RZ=A33 
PEAR=0 
CE>AR  =  0 
R6AR=C 
APGUN= 
G  A  f  M  A  R 
GANMA= 
VXE=A1 
VVB=A2 
VZB=A3 
VX1=A1 
VY1=A1 
VZl=Al 

vzxe=v 

VYVB=V 
ALPhAR 
BETAR= 
ALPhAD 
BETAD= 
S ALPHA 
CALPFA 
SL;ETA  = 
CEETA= 

R  E  T  U  R  l\ 

ENC 


TINE 

IT  R 

/SFI 

LTAR 

AR,X 

A23  , 

i  RES 

ALPH 

R,XF 

f* ,  /.  u 

=  THE 

HI/5 

SI/5 

=  CSI 

=  DCC 

SIN( 

CCS( 

SIN{ 

CCS( 

TFET 

HI*C 

Hi-; 

PET 

HI* 

H  I* 

HETA 

Hl-S 

HI*S 

.5*R 

AMBO 

T(VS 

GAMD 

GAMB 

-•-GAM 

.5*P 

.5*C 

.  5*R 

GAMC 

=  BAP 

GAMM 

1*;GA 

1*GA 

1*GA 

i*VX 

2*VX 

3*VX 


zb/v 

YB/V 
=  DAT 
OARS 
=  ALP 
BETA 
=  CSI 
=  DCC 
DSIN 
CCCS 


CONST 
EAL*3(A 
N/CAMUG 
t THRUST 
I  ,  Y  I  ,  l  I 
VS«J,VX1 
(to)  ,ALP 
A,SGAivM 
CRCE.XF 
CT, NPRI 
TA/57.2 
7.29577 
7.29577 
N(THETA 
S(THETA 
PHIR) 
PHIR) 
PSIR) 
FSIR) 
A 

THETA 
THETA 
*  S  P  S I 
THETA*S 
THETA* S 
*CPSI 
THETA-C 
THETA*C 
HO*S 
T*R)*(G 
G) 
GT 
CT 
OCT 
*B/V 
*CBAR/V 
Z-B/V 
CT*R/V 
CCStARG 
AR*57.2 
MDOT-R+ 
MDOT*R+ 
KDGT-P.+ 
b+A21*V 
B+A22*V 
E+A23-V 
XB 


psi+cpfi*cpsi 
psi-sphi*cpsi 

psi-cphi*spsi 
psi+sphi-spsi 

amcgt*r)+zcct*zcct+rdot*rcct 


UM) 

55779 

A12*RD0T+A13*ZCCT 

A22*RCCT+A23*ZCCT 
A32*RD0T+A33*ZCCT 

YE  +  A3r-VZB 
YB+A32-VZB 
YB+A33-VZB 


AN( VZX3) 
IN(VYVB) 
HAR*57. 295779 
R-57. 255779 
N( ALPHAR) 
S( ALPhAR) 
(  b  E  T  A  R ) 
(bETAh  ) 
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SLBPCLTINE  CLCKUP 


PIRPCSE 


C** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


TAKES  THE  PASSED  VALUES  CF  ANGLE-CF 
AND  SIDESLIP  AND  PERFORMS  A  TtoO-hAY 
INTERPOLATION  CF  ALL  lb     CCEFAB  C  /i  L  P 
BETA  DEPENDENT  COEFFICIENTS  )  ARRAY 
ATIVBS.   A  LINEAR  INTERPOLATION  IS 
CN  THE  9  CCEFA  (  ALFFA  CNLY  DEPENCE 
FICIENTS  )  ARRAY  DERIVATIVES.   TFE 
LATEC  COEFFICIENT  VALUES  (  CR(I).I= 
STORED  IN  COM  NGN  MEMORY  FCP  USE  BV 
SSEQMS  IN  EVALUATING  THE  EQUATIONS 


■if 

-ATTACK  * 
LINEAR  * 
FA  AND  * 
CERIV-  * 
PERFORMED  * 
NT  COEF-  * 
INTERPG-  * 
1,24  )  ARE* 
SUBROUTINE* 
CF  MGTICN.* 


:**s{5!<*********** 


c 

c 


SL 

IN 

CO 

IT  A 

2S, 

21, 

4CF 

5LP 

6,C 

7C  N 

AI 

El 

IA 

ie 

ia 

IE 
DX 
CY 


EROU 
PLIC 
MMGN 
E,OE 
B,Cb 
A32  , 
(24) 
FArC 
PHI, 
,ZMC 
=  0.2 
=  0.1 
1=AI 

1  =  31 

2  =  IA 
2=IB 
-bl- 
=  AI- 


12,A13,A21  ,A22,A23, A3 
FLJAR,UBAR,RBAR,FRhCS, 


TINE    CLCKUP 

IT    REAL*8(A-h,G-Z) 

/SPIN/GAhDGT,GAMMA,RDCT,PFI  , THETA,PSI  ,  DELTA A,  DEL 

LTaR,THF.UST,V,C0EFAB(15,2G,1Q},C0EFA<2C,9)  , SSAM, 

ARtXI,Yi,ZI,ZXI,G,RHC»All,A 

A23,VSQtVXl,VYl,VZlf P,G,RZ, 

,RES<6),ALPPACtALPHAR,BETAC,BETAR,SEFTA,CBETA,SA 

ALPHA,SGaMMA,CGAMNA,SPSI ,CFSI,STHETA,CTPiETA,SPPI 

R,XFGRCE,XFCR1, Y FORCE fY FORI, ZFORCE,Z FORI, XMOM,YM 

y,ZDCT,NPRINT,NANE(4) 

*ALPHAD+1.G 

*BETAD+5.0 


1  +  1 
1  +  1 

iei 

IA1 


INTERPOLATE  ANGLE-GF-ATT ACK 
DEPENDENT  COEFFICIENTS 


ANC  SIDE SL  IF 


CO  1  K=l,  16 

C1=C0EFAB(K, IA1,  IB  1  )  +  CX* I  CO EFAB ( K ,  I A  1,  IB2  )- 
lCCEFABi  K,  IA1,  IblJ  ) 

C2=C0EFAb(K,IA2,  IB  1) +  DX*< COEFAB ( K ,  IA2fIB2)- 
1CCEFAB( K,  IA2,  IB1  )  ) 

1  CR(K)=C1+DY-(C2-Clf 

INTERPOLATE  ANGL E-OF-ATT ACK  CEFENDENT  CECFFICIENTS 

CO  2  M=l,9  • 
I=N+15 

2  CR(  I)=CC;EFA(  IAi,  M)+DY*  (COEF  A  ( I  A 2  ,  P  )  -COE F  A  (  I  A  1 , M  )  J 
RETURN 

ENC 
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PLP 


Q  4  $  4  4  $ 

C   SLB 
C 

C 

c 

c 
c 
c 
c 

c 


PCUTIN 
PCSE 


E  SSECNS 


:****j(>*******  ******$:$:*$;»:**  **£****# 


CALC 
F I C  I 
CF  T 

ARE 

steR 

FUNC 


ULATE 
EMS 
HE  SI 

MADE 
CUT  IN 
TION. 


S  THE  TOTAL  FORCE  AND  MCNENT  COEF 
ANC  THE  RESIDUALS  (  RES(I),I=1,6 
X  EQLATIONS  OF  NOTION.   TEE  RESID 
AVAILABLE  IN  CCNNCN  MEMCPV  FOR 
E  FN  TO  EVALUATE  TFE  CRITERIA 


A  *  ** 
* 
* 
* 

)    * 

UALS 


* 
* 

Q*^*^>*?5  **-;***:?*******************  **J?***^******^^*^***  A*  ***^^ 


C 
C 
C 


C 
C 

c 


SLBROUTINE  SSEQNS 

IMPLICIT  REAL*3( A-F,0-Z) 

CCNNCN/SPIN/GAMDCT,GAN,NA,RDCT,PFI  ,TPETA,PSI,DELTAA,DEL 
1TAE,DELTAR,THRUST,V,COEFAB(15,20,1C),COEFA(2C,9),SSAN, 
2S,B,CBAR,XI,YI,ZI, ZXI , G ,RhO , Al 1 , A  12 , A13 , A 21, A  22, A23,A3 
31,A32,A33,VSQ,VXl,VYl,VZl,P,C,RZ,FEAR,J6AR,FdAR,FRhCS, 
4CP(24 ) ,RES(6 )  f ALP HAD, ALPHAk , BET  A C , BET AR , SB E TA , CBE T A  ,  SA 
5LPhA,CALPHAfSGANNA,CGANNA,SFSi,CPSI,STHETA,CTFETA,SPFI 
6,CPHI,R,XF0RCF,XFCRL, Y FORCE, Y FORI, I  FORCE , ZFCR 1 , XNCN , YM 
7GN,ZMCN,ZDCT,NPRINT,NAME(4) 

CALCULATE  THE  TOTAL  FORCE  ANC  NONENT  CCEFFICIEN7S 


*DELTAE+CR(24)*QBAR 

L6)*PBAR+CR(  19P>RBAR  +  CR<  10J -CELTAA  + 


CX=CR(5  )+CR 
CY=CR(4)+CR 
ICR(7)*DELTAR 
CZ=CR(6)+CR(13)*DELTAE*CR(22)*QEAR 

CL=CR(1 )+CR(ll)*DELTAA*CR(8)*0ELTAR+CR( 17)*P8AR+ 
1CR(20)*RBAR 

CN=CR (2  )+CR(  14 ) *D ELT AE+CR{ 23 )*Q BAR 

CN  =  CR(3)+CR(  12)*CELTAA+CR(9)*DELTAR+CR(  18MPBAR  + 
lCR.(21)*RbAR 

C=1.0-ZXI*ZXI/(XI*ZI ) 


X  FORCE 
YFCRCE 
ZFCRCE 
XNCN^H 
YMCM=H 
ZNCN=H 
XFCR1= 
YFORl= 
2FCRL= 


CALCULATE  THE  AERCCYNAMIC  FCRCES  AND  NCNENTS 

+TFROST 


=  PRHC 
=  hRHO 
=  PRHO 
RFCS* 
RhOS  ': 
RFCS* 
A  L 1  *  X 
A  12*X 
A13*X 


S* 
S* 
S  * 
VS 

vs 
vs 

EC) 

FO 
FO 


VSQ 
VSu 
VSQ 

e-B 

Q*C 
Q*B 

KCE 
RCE 
RCL 


*cx 

*CY 

:!:CZ 
*CL 
BAR 
*CN 
+A2 
+  A2 
+  A2 


*CM 

1*YFCRCE+A31*ZFCRCE 

2*YFCRCE+A32*ZFCRCE 
3*YFORCE+A33iZFCRCE+SSAN*G 


RES 
RES 
PES 
RES 

1  )/Z 

2/C 
RES 

1  (PZ 
PES 

11.) 
RET 
END 


CALCULATE  THE  EQUATION  RESIDLALS 

(1)=>F0R1/<SSAN^V) 

<2)=<YFGR1-SSAN.*GANCCT*GAMD0T*R)/<SSAN*V  ) 
(3)=ZF0R1/(SSAN*V) 

<4)=(((YI-ZI)/XI-(ZXI*ZXI)/{XI*ZI))*Q*PZ-M1.-(YI-XI 
I  )*(  ZXI*P*Q/XI  )-KHRFOS*VSQ*B/XI  )  * (CL+ < ZX  I*CN/ZI))  ) 

(5)=HRHOS*VSQ*CBAR*CM/YI+(ZI-XI )*P*RZ/Y I-M  ZXI/ YI )* 
*RZ-F*P ) 

<c>  =  UZXI*ZXI/(XI*ZI)-(YI-XI)/ZI)*P*C+UYI-ZI)/XI- 
*ZXI*Q*RZ/ZI+<HRHOS*VSQ*B/ZI )*(ZXI*CL/XI+CM)/D 

URN 
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F,   SAMPLE  OUTPUT 


CONFIGURATION  A 

-  INITIAL  CCNDITICNS  - 

GANCCT  R  RDCT  PHI         TPETA 

(RAC/SEC)  (FT)  (FT/SEC)  (DEG)        (LEG) 

1.01000  23.CC000  0.0  5.00000   -35.00CCC 

PSI  CELTAA  DELTAR  DELTAE      TFBUST 

(CEG)  (DEG)  (DEG)  (DEG)       (LBS) 

75.CC000  15.00COO  -30.00000  -25.00000     0.0 

ZDCT  RFC  DFMX  OXMAX       FMAX 
(FT/SEC) 

340.G0C00  C.C0089  0.1000D-60  C.ICOOC-60  0.1000D-6C 

LNAX  IV 

-1C5C  2 


-  EXTREN  ITERATICNS  -  -  - 

STAGE  NO.  ttt  TRIAL  NO.  11830 

DL=     0.0 

FUNCTION    \/4LUE=  0.25726C-23 


AP( 

1)  = 

C  .S8C60D 

OC 

DSC 

1)  = 

0.222040-16 

AP  ( 

2)  = 

C  .698780 

02 

DS( 

2)  = 

0.542  1CD-21 

AR( 

3)  = 

C. 343450 

C3 

DS( 

3)  = 

C.33881D-19 

AR( 

A)  = 

0  .202990 

02 

0S( 

A)  = 

C.  14*110-15 

AR( 

5)  = 

-C  .39  L63D 

02 

0S( 

5)  = 

0.  1  L^b9D-i4 

AR( 

6)  = 

0.217760 

C2 

DS( 

£)  = 

0.847030-20 
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CGNF  IGLRATICN  A 


-  -  -  FINAL  PARAMETERS  -  -  - 


ALPHAD 
(DEG) 

4£. 38819 


PS  I 
(DEG) 

69.87808 


EETAC 
(CEG) 

1C  .457L9 


GAMMA 
(CEG) 

£4  .81609 


GAMCOT 
(RAD/SEC) 


ThETA 
(DEG) 


0.98060   -39.16334 


R 
(FT) 

21.77595 


THRUST 
(LBS) 

0.0 


PFI 
(DEG) 

20.29875 
RHO 
C.CCC£9 


(FT/SEC) 
344.86437 


VX1 
(FT/SEC) 

31.15947 


VY1 

(FT/SEC) 


VZ1 
(FT/SEC) 


•C.COCCC   343.45381 


(RAD/SEC) 
0.61928 


Q  PZ 

(RAD/SEC)  (PAC/SEC) 


C. 26376 


C. 71309 


RCCT 
(FT/SEC) 

0  •  c 


XFORCE= 

VFORCE= 

2FCRCE= 

XMOM  = 

YNCM= 

ZMCM  = 


2992.653C0 

-7896.07153 

-68433.22273 

5427.36427 

-127703 .55  12  1 

42513.27919 


RES(1)= 

PES(2)= 
RES(2)= 
RES(4)= 
FES(5)= 
RES(6)= 


-0.760C5D-12 


-0.13141C-12 


0. 13899D-11 


-0.32226D-13 


0.38337C-13 


0.  1U54C-12 


SIM  CF  ThE  Su LAKES 
CF  THE  RESICLALS 


0.25726C-23 
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G.   PROGRAMMING  CONSIDERATIONS 

Program  SPIN  requires  66K  bytes  of  storage.   When  buffer 
space  is  taken  into  account,  a  request  for  80K  bytes  should 
be  made  when  submitting  the  job  for  processing. 

The  time  requirement  is  variable  and  is  very  dependent 
on  the  closeness  of  the  initial  conditions  to  the  steady 
spin  condition.   Convergence  has  been  obtained  with  a  total 
CPU  time  of  a  little  over  one  minute  to  as  high  as  15  min- 
utes; this  time  includes  about  25  seconds  for  the  compile 
and  link  steps. 

The  optimization  scheme  loop  time  is  very  rapid.   The 
time  required  to  complete  one  stage  is  0.2  second,  which 
requires  18  criteria  function  evaluations.   The  criteria 
function  evaluation  occurs  at  the  rate  of  one  per  .0111 
second  and  this  accomplishes,  as  Figure  Bl  indicates,  the 
evaluation  of  the  equation  constants,  interpolation  of  the 
24  aerodynamic  coefficients,  and  the  evaluation  of  the 
equations  of  motion.   By  looking  at  the  listing  for  each 
of  the  subroutines  involved,  one  can  better  appreciate  the 
speed  with  which  the  program  operates. 

The  variable  time  requirement  of  the  program  can  create 
a  problem  to  the  researcher.   More  often  than  not,  the 
turn  around  time  for  any  program  is  determined  by  the  time 
and  storage  requested.   As  a  general  guideline  to  the  use 
of  program  SPIN,  the  following  is  suggested: 

a.   Set  LMAX  (the  maximum  number  of  stages)  to  a  value 
of  1,050. 
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b.  Request  80K  storage  and  a  four  minute  run  time. 

c.  Check  results  of  the  run.   If  the  value  of  the 

—5 
criteria  function  is  greater  than  10   ,  guess 

another  set  of  initial  conditions  and  rerun.   Re- 
sults indicate  that  if  the  value  of  the  criteria 
function  hasn't  dropped  below  the  above  value  in 
the  four  minutes  alloted  that  an  increase  of  time 
will  not  yield  a  steady  spin  condition.   The  pro- 
gram in  this  case  is  converging  on  a  local  mini- 
mum which  is  not  indicative  of  a  steady  spin.   If 

the  value  of  the  criteria  function  is  less  than 

—5 
10   ,  rerun  the  program  with  LMAX  set  at  2,500 

and  a  requested  time  of  nine  minutes  if  greater 

accuracy  is  desired. 
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APPENDIX  C 

AIRCRAFT  DATA 

The  following  aircraft  data  was  used  in  the  study  of 
steady  spin  conditions.   Table  CI  tabulates  the  mass  and 
dimensional  characteristics  of  the  two  aircraft  and  the 
following  pages  list  the  aerodynamic  coefficients  for  each 
configuration . 

TABLE  CI.   Mass  and  Dimensional  Characteristics 


Characteristics 

Configuration  A 

Configuration  B 

m  (slugs) 

1554.00 

771.81 

s  (ft2) 

525.00 

695.00 

b  (ft) 

63.00 

38.00 

c  (ft) 

9.04 

23.76 

2 

I  (slugs— ft  ) 

X 

53100.00 

13600.00 

2 
I  (slugs— ft  ) 

y 

299000.00 

128000.00 

I  (slugs— ft  ) 
z     r 

338750.00 

138000.00 

I   (slugs-ft2) 
xz     r 

12480.00 

4340.00 
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